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Introduction 


The  meshless  method  implemented  in  PHLEXcrack^^  is  the  so-called  Element  Partition 
Method  (EPM)  which  is  a  variation  of  /ip-Cloud  [15-17, 30,  38]  and  Finite  Element  Parti¬ 
tion  of  Unity  [2,3,33]  Methods.  The  key  feature  of  the  Element  Partition  Method  is  the  use  of 
a  partition  of  unity  to  build  the  approximation  spaces.  This  partition  of  unity  framework  has 
several  powerful  properties  such  as  the  ability  to  produce  seamless  hp  approximations  with 
nonuniform  h  and  p,  the  ability  to  develop  customized  approximations  for  specific  applications 
like  dynamic  crack  propagation,  the  capability  to  build  p-orthotropic  approximations  on,  e.g., 
three-dimensional  tetrahedral  meshes,  etc.  Several  of  the  so-called  meshless  methods  proposed 
in  the  last  years  also  make  use  explicitly  or  implicitly  of  a  partition  of  unity  to  build  the  ap¬ 
proximation  spaces  [4,17,31,35].  The  fundamental  difference  between  these  methods  and  the 
EPM  is  in  the  choice  of  the  partition  of  unity.  In  the  EPM  the  partition  of  unity  is  provided 
by  a  combination  of  Shepard  [25,28,45]  and  finite  element  partitions  of  unity.  This  partition 
of  unity  allows  the  inclusion  of  arbitrary  crack  geometry  in  a  model  without  any  modification 
of  the  initial  discretization.  We  call  this  partition  of  unity  Cracked  Element  Partition  of  Unity 
(CEPOU).  This  choice  of  partition  of  unity  also  avoids  the  problem  of  integration  associated 
with  the  use  of  moving  least  square  or  conventional  Shepard  partitions  of  unity  which  are 
used  in  several  meshless  methods  [4,9,27,29,31,35].  Performance  studies  show  that  the  EPM 
is  orders  of  magnitude  more  computationally  efficient  than  other  meshless  methods.  In  fact, 
the  same  studies  also  show  that  the  EPM  is  often  more  computationally  efficient  than  the 
Finite  Element  Method,  while  being  flexible  enough  to  allow  the  modeling  of  dynamic  crack 
propagation  in  complex  three  dimensional  structures  without  any  remeshing.  In  addition  to 
that,  the  partition  of  unity  used  in  the  EPM  allows  ease  implementation  of  essential  boundary 
conditions  which  is  a  problem  when,  e.g.,  moving  least  square  partitions  of  unity  are  used  as 
in  the  Element  Free  Galerkin  and  Reproducing  Kernel  Methods  [4,31,35]. 

An  important  feature  of  the  EPM  is  the  ease  in  which  customized  shape  functions 
can  be  added  to  the  set  of  interpolating  functions.  This  feature  can  be  very  advantageous  in 
cases  where  information  about  the  solution  and  how  it  looks  is  known  a  priori,  as  in  the  case 
of  re-entrant  corners.  For  cracks,  asymptotic  singular  functions  are  used  to  enrich  the  basis 
functions  of  the  meshless  nodes  close  to  the  crack  tip  with  the  order  of  singularity  derived 
from  analytical  solutions  [13, 14].  Using  customized  shape  functions  reduces  the  number  of 
degrees  of  freedom  needed  to  accurately  represent  the  solution  thus  reducing  time  and  storage 
requirements. 

Once  a  solution  is  obtained  at  some  step,  the  amount  and  direction  of  crack  propagation 
over  the  next  time  increment  can  be  estimated.  These  values  are  evaluated  at  a  series  of 
representative  points  along  the  crack  front.  The  front  is  then  advanced  to  its  new  position. 
The  model  is  updated  to  accommodate  the  advanced  crack.  However,  unlike  with  finite 
elements  the  crack  can  propagate  freely  in  any  direction  and  not  necessarily  along  element 
faces.  Moreover,  if  the  crack  cuts  an  integration  cell  it  is  not  required  to  regenerate  any  new 
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cells  or  connectivities.  However,  for  an  accurate  representation  the  integration  rule  will  be 
enriched/modified  within  those  cells  that  are  cut  by  the  crack.  In  addition,  enrichment  with 
the  auxiliary  singular  functions  will  be  maintained  for  the  meshless  nodes  closest  to  the  crack 
front. 

The  crack  propagation  problem  usually  requires  complex  geometric  and  topological 
descriptions  and  calculations.  This  serves  several  essential  tasks  such  as.  defining  the  crack 
surface  and  front,  updating  them,  interaction  of  crack  with  outside  boundary,  discontinu¬ 
ities  and  integration  around  crack,  etc.  In  order  to  obtain  accurate  and  efficient  geometric 
calculations,  the  geometric  kernel  of  HyperMesh^^  (a  powerful  geometry  engine  and  pre/post¬ 
processor)  was  modified  for  our  purposes,  and  appropriate  interface  functions  were  designed 
and  coded.  This  kernel  is  being  linked  as  a  library  into  PHLEXcrack  .  PHLEXcrack  code 
is  based  on  a  proprietary  FEM  environment  ProPHLEX,  which  was  during  this  project, 
expanded  to  handle  the  EPM  and  CEPM  methods  [19]. 


2  Formulation  of  Element  Partition  Approximations 

Let  Q.  be  an  open  domain  in  iE”,  n  =  1, 2, 3  and  Tn  an  open  covering  of  Q  consisting  of  N 
supports  uja  (often  called  clouds)  with  centers  at  x  Q;,  l,...,iV,  i.e., 

N 

Tn  =  C  [J  Wq, 

a=l 


The  basic  building  block  of  an  element  partition  of  unity  (EPOU)  approximation  is 
a  set  of  functions  (pa  defined  on  the  supports  a  =  1,...,IV,  and  having  the  following 
property 

Pa  e  CIM,  s  >  0,  1  <  a  <  iV 

Y^(Pc,{x)  =  l 
a 

The  last  property  implies  that  the  functions  Pa,  a  =  1, ...  ,N ,  axe  non-zero  only  over  the 
supports  Wq,  O'  =  1, . . . ,  A^.  The  functions  Pa  are  called  a  partition  of  unity  subordinate  to 
the  open  covering  Tn-  Examples  of  partitions  of  unity  are  Lagrangian  finite  elements,  moving 
least  squares  and  Shepard  functions  [17,25]. 


2.1  Finite  Element  Partition  of  Unity 

In  the  case  of  finite  element  partitions  of  unity,  the  supports  (clouds)  are  simply  the  union 
of  the  finite  elements  sharing  a  vertex  node  Xa  (see,  for  example,  [33, 38]  and  Figs.  1  and  2) 
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and  the  partition  of  unity  function  ipa  is  equal  to  the  global  finite  element  shape  function  Na 
associated  with  a  vertex  node  as  depicted  in  Fig.  1.  In  this  case,  the  implementation  of  the 
method  is  essentially  the  same  as  in  standard  finite  element  codes,  the  main  difference  being  the 
definition  of  the  shape  functions  as  explained  in  Section  2.3.  This  choice  of  partition  of  unity 
avoids  the  problem  of  integration  associated  with  the  use  of  moving  least  square  or  Shepard 
partitions  of  unity  used  in  several  meshless  methods.  Here,  the  integrations  are  performed 
with  the  aid  of  the  so  called  master  elements,  as  in  classical  finite  elements.  Therefore,  the 
EPM  can  use  existing  infrastructure  and  algorithms  for  the  classical  finite  element  method. 


Figure  1:  Global  finite  element  shape  function  built  on  a  mesh  of  triangles  and  quadrilat¬ 
erals. 


2.2  Shepard  Partition  of  unity 

In  this  section  we  describe  the  construction  of  a  partition  of  unity  using  the  so  called  Shepard 
formula. 

Let  Wa  denote  a  weighting  function  with  compact  support  Ua  that  belongs 

to  the  space  CQ{oJa),  s  >0  and  suppose  that 

Wa(x)  >0  V  X  G 

Suppose  that  such  a  weighting  function  is  built  at  every  cloud  Ua,  a  =  l,...,N.  Then 
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Figure  2:  Overlapping  patches  corresponding  to  clouds  u)a  and  Polynomials  of  differing 
degree  Pa  and  can  be  associated  with  nodes  at  Xa  and  Xp  so  as  to  produce  non-uniform 
hp-cloud/finite  element  approximations. 


the  partition  of  unity  functions  (pa  associated  with  the  cloud  u)a  can  be  defined  by 


f  \  - 


^  6  {7  :  >V-y(x)  ^  0} 


(1) 


which  are  known  as  Shepard  functions  [45]. 

Consider  now  the  case  in  which  the  weighting  function  Wq  is  taken  as  the  global  finite 
element  shape  functions  Na  associated  with  the  node  Xa,  a  1,. .  .,N.  Let  r  be  a  finite 
element  with  nodes  x^,  /5  G  Jr  where  Ir  is  an  index  set  and  let  x  €  r .  The  only  non-zero 
partition  of  unity  functions  at  x  are  given  by 


(Pa{x) 


Nqjx) 

E/3eIr 


Ngjx) 

1 


Ng{x) 


CX  ^  Xr 


sincG  the  finite  element  shape  functions  form  a  partition  of  unity.  Therefore,  it  is  not  necessary 
to  use  the  Shepard  formula  to  build  the  partition  of  unity  when  the  weighting  functions  are 
taken  as  global  finite  element  shape  functions.  However,  as  demonstrated  below,  the  above 
formula  is  the  key  to  build  a  partition  of  unity  when  the  finite  element  r  is  severed  by  a  crack. 
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2.3  Element  Partition  Shape  Functions:  The  Family 

The  construction  of  the  so-called  element  partition  shape  functions  is  based  on  the  following 
observation: 

Let  {Li}i^x  denote  a  set  of  functions  which  can  approximate  well,  in  an  appropriate 
norm  the  solution  u  of  a  boundary  value  problem,  i.e.,  there  exist  Ui,  i  £l  such  that 


where  Uhp  =  Eiei  ^  denotes  an  index  set. 

Now  consider  the  following  set  of  cloud  or  element  partition  shape  functions,  defined  as 

0“  :=  (paLi,  a  =  l,...,N,  i  el 

where  cpa  is  a  partition  of  unity  function.  Then  it  is  not  difficult  to  show  that  linear  combina¬ 
tions  of  these  EP  shape  functions  can  also  approximate  well  the  function  u 


ot  i 


51/  Ui^aLi  —  ^  ^  UiLi 

a  i  a  i 

^  5^  ^a'^hp  ^  '^hp  5  V  “  '^hp 

a  OL 


(2) 


Note  that: 

(i)  The  EP  shape  functions  ,  a  =  1, . . . ,  iV,  i  €  X,  are  non-zero  only  over  the  cloud  Wo,. 

(ii)  Linear  combination  of  the  EP  shape  functions  can  reproduce  the  approximation  Uhp  to 
the  function  u. 

(iii)  The  functions  Lj,  *  €  X,  can  be  chosen  with  great  freedom.  The  most  straightforward 
choice  is  polynomial  functions  since  they  can  approximate  well  smooth  functions.  How¬ 
ever,  for  many  classes  of  problems,  including  the  case  of  fracture  mechanics  problems, 
there  are  better  choices.  This  case  is  discussed  in  details  in  Section  3. 

In  this  section,  element  partition  shape  functions  are  defined  in  an  n-dimensional  setting 
using  the  same  ideas  outlined  above. 

Let  the  functions  a  =  1, . . . ,  iV,  denote  a  finite  element  partition  of  unity  subor¬ 
dinate  to  the  open  covering  7iv  =  {^0)0=1  ^  domain  Q,  c  iK",  n  =  1, 2,3.  Here,  N  is  the 

number  of  vertex  nodes  in  the  finite  element  mesh.  The  cloud  tUa  is  the  union  of  the  finite 
elements  sharing  the  vertex  node  Xa- 
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Let  Xa(wa)  =  span{LiQ}igx(a)  denote  local  spaces  defined  on  a  =  where 

I  {a),  a  =  l,...,N,  axe  index  sets. 

Suppose  that  the  finite  element  shape  functions,  (fa,  are  linear  functions  and  that 
Vp-iM  C  a  =  l,...,N, 

where  Vp-i  denotes  the  space  of  polynomials  of  degree  less  or  equal  to  p  -  1.  The  element 
partition  shape  functions  of  degree  p  are  defined  by  [33, 38] 

=  {<pf  =  ^aLia,  a  =  1, . . . ,  iV,  i  e  I{a)}  (3) 

Note  that  there  is  considerable  freedom  in  the  choice  of  the  local  spaces  Xa-  The 
most  obvious  choice  for  a  basis  of  Xa  is  polynomial  functions  which  can  approximate  well 
smooth  functions.  In  this  case,  the  EPM  is  essentially  identical  to  the  classical  FEM.  The 
implementation  of  hp  adaptivity  is,  however,  greatly  simplified  by  the  EP  framework.  Since 
each  basis  functions,  {Liah^x^a),  a  =  1, . . . ,  AT,  can  have  a  different  polynomial  order,  we 
can  have  different  polynomial  orders  associated  with  each  vertex  node  of  the  finite  element 
mesh  [38].  The  approximations  can  also  be  non-isotropic  (i.e.,  different  polynomial  orders  in 
different  directions),  regardless  of  the  choice  of  finite  element  partition  of  unity  (hexahedral, 
tetrahedral,  etc.).  The  concept  of  edge  and  middle  nodes,  which  is  used  in  conventional  p 
FEMs,  is  not  needed  in  the  framework  of  EPM.  The  implementation  of  h  adaptivity  is  also 
simplified  in  EP  methods  since  it  needs  to  be  done  only  on  the  partition  of  unity  (linear 
finite  elements  in  this  case).  Therefore,  the  implementation  of  h  adaptivity  for  high-order 
approximations  is  the  same  as  for  linear  approximations  (there  is  no  need,  for  example,  of 
using  high-order  constraints  as  done  in  traditional  hp  finite  element  methods). 

There  are  many  situations  in  which  the  solution  of  a  boundary  value  problem  is  not 
a  smooth  function.  In  these  situations,  the  use  of  polynomials  to  build  the  approximation 
space,  as  in  the  FEM,  may  be  far  from  optimal  and  may  lead  to  poor  approximations  of  the 
solution  u  unless  carefully  designed  meshes  are  used.  In  the  EPM,  we  can  use  any  a-priori 
knowledge  about  the  solution  to  make  better  choices  for  the  local  spaces  Xa-  For  example,  in 
Section  13,  we  solve  a  boundary  value  problem  in  which  the  solution  possesses  point  or  lines 
of  singularities  at  some  parts  of  the  domain.  Then  we  can  use  the  local  spaces  Xa  to  build 
element  partition  shape  functions  that  represent  these  singularities  much  more  effectively  than 
polynomial  functions.  The  construction  of  these  customized  shape  functions  is  discussed  in 
details  in  Section  3. 

Suppose  that  the  partitions  of  f2  into  finite  elements  which  form  the  partition  of  unity 
satisfy  the  usual  regularity  assumptions  in  two  and  three  dimensions,  and,  in  one  dimension, 
the  ratio  between  the  lengths  of  neighboring  elements  is  bounded.  We  denote 

ha  =  diam(a;a) 

h  =  max  ha 

a=l,...,N 
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and 


=  span{^f },  cx  =  l,...,N,  i  e  X{a) 

where  the  EP  shape  functions  0“  are  defined  in  (3).  In  addition,  suppose  the  unity  function, 
1,  belongs  to  the  set  of  local  approximation  functions,  i.e., 

1  £  XQ!(^a)  ^  1>  •  •  •  5 


and  that  there  exists  a  quantity  e  depending  of  a,  h,  p,  and  u,  such  that 

||u  —  ^  l,...,iV 

Then,  it  can  be  proved  that  [17, 33]  there  3  Uhp  6  X^^  such  that 


/N(h)  ^ 

V  Q=1  J 


1/2 


where  the  constant  C  is  independent  of  u,  h,p. 


2.4  Construction  of  Discontinuous  Element  Partition  Shape  func¬ 
tions 

In  this  section,  a  technique  to  modify  a  finite  element  POU  at  elements  cut  by  a  crack  surface  is 
described.  The  POU  is  modified  such  that  a  discontinuity  in  the  displacement  field  across  the 
crack  surface  is  created.  The  POU  is  modified  only  for  elements  cut  by  the  crack.  Elsewhere, 
a  finite  element  partition  of  unity,  as  described  in  Section  2.1,  is  used.  The  technique  allows 
for  elements  to  be  arbitrarily  cut  by  the  crack  surface  without  any  mesh  modification.  In 
these  elements,  the  POU  is  built  using  the  Shepard  formula  (1)  and  finite  element  shape 
functions  as  weighting  functions.  The  technique  is  described  in  a  two-dimensional  setting  in 
order  to  simplify  the  notation  and  pictures,  however,  the  technique  can  be  used  without  any 
modifications  in  the  three  dimensional  case.  Three-dimensional  problems  are  solved  in  Section 
13  using  this  technique.  We  also  restrict  ourselves  to  the  case  of  quadrilateral  elements.  The 
case  of  triangular  elements  is  identical. 

The  technique  is  described  using  the  discretization  depicted  in  Fig.  3.  For  elements 
that  do  not  intersect  the  crack  surface,  the  POU  is  provided  by  the  finite  element  shape 
functions.  For  example,  at  any  point  x  G  n,  there  are  four  partition  of  unity  functions 

(fi  =  Ni,  (p2  =  N2,  (p5  =  Ns,  (pe  =  Ns 

Since  the  FE  shape  function  form  a  POU,  the  following  holds 

+  (p2{x)  +  <P5(®)  +  ¥’6(®)  =  Ni{x)  -t-  N2{x)  -f-  iV5(x)  +  Ns{x)  =  1  V  X  G  ri 
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Figure  3:  Example  of  element  partition  discretization  cut  by  a  crack 
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The  Shepard  formula  (1)  can  also  be  used  to  build  the  POU  but  in  the  case  of,  for  example, 
element  n,  the  results,  in  this  case,  are  trivial,  i.e.. 


(Pa{x) 


_ _ 

Ni{x)  +  N2{x)  +  TVs  (a;)  +  N^ix) 


Naix) 


O'  =  1,2, 5, 6 


For  the  elements  that  are  cut  by  the  crack  surface,  the  Shepard  formula  (1)  is  used  in  combi¬ 
nation  with  the  visibility  approach  [4].  That  is,  the  crack  surface  is  taken  as  an  opaque  surface 
and,  at  a  given  point  x  er,  the  weighting  function  Na  is  taken  as  non-zero  if  and  only  if  the 
segment  connecting  x  and  Xa  does  not  intersect  the  crack  surface.  As  an  illustration,  consider 
the  case  of  the  element  T2  shown  in  Fig.  3.  At  the  point  y  e  T2,  according  to  the  visibility 
criteria,  the  only  non-zero  weighting  functions  are  Ne  and  Nr,  since  the  segments  [y  -  Xw] 
and  [y  -  ®ii]  intersect  the  crack  surface.  The  partition  of  unity  at  y  €  r2  is  then  given  by 

/  ^  _  ^6(2/)  ^  Nr{y) 

“  NM  +  N,{y)  Ne(y)  +  N,(y) 


Note  that 


^6{y)  +  My)  = 


Nsiy)  +  Nriy) 


=  1 


Nsiy)  +  Nr{y) 

Therefore,  the  functions  (pe  and  (pr,  as  defined  above,  constitute  a  partition  of  unity. 

At  the  point  z  e  T2,  the  partition  of  unity  is  given  by 

^  ^  _  Ario(;2;)  ^  Nnjz) 

+  ^  N,o{z)  +  Nn{z) 


Therefore,  the  approximation  is  discontinuous  across  the  crack  surface. 

As  another  example,  consider  the  case  of  the  element  T3  with  nodes  X5,  *6,  and  Xiq 
as  depicted  in  Fig.  3.  At  the  point  r  G  ra,  the  partition  of  unity  is  given  by 


TV5(r) 

iV5(r) -t- iV6(r)  +  iV9(r) 


Mr)  = 


iVeCr) 

iV5(r)  +  iV6(r)-hiV9(r) 


'  N,(T)  +  Ni{r)  +  N,(T) 

At  the  point  s  G  ra,  the  partition  of  unity  is  given  by 


(^io(s) 


TVio(s) 

iVio(s) 


1 


Let  us  now  show  that  the  approximation  is  continuous  at  points  not  at  the  crack  surface. 
Consider,  for  example,  the  point  t  located  at  the  boundary  between  elements  T2  and  tq,  as 
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M 


\ 


shown  in  Fig.  3.  Consider  the  POU  function  (pio  associated  with  node  xio-  If  this  function  is 
computed  from  element  T2  we  have 


iVio(t) 

Niait)  +  Niiit) 


If  the  function  ipiQ  is  computed  from  element  re  we  have 


_ Nw{t) _  Niojt) 

Nio{t)  +  +  Nis{t)  +  Nx4{t)  Nio{t)  +  Nu{t) 


Therefore  ^io^(^)  =  ^io^(^)  and  the  function  (pio  is  continuous  at  t.  This  same  argument  can 
be  used  at  other  points  not  located  at  the  crack  surface. 

At  the  elements  intersecting  the  crack  front,  like  element  T4  shown  in  Fig.  3,  the  crack 
is  modeled  with  the  aid  of  singular  functions  as  described  in  Section  3. 

Note  that  the  above  technique  allows  arbitrary  cut  of  the  finite  element  mesh  by  the 
crack  surface.  Also,  the  computational  cost  of  EP  shape  functions  is  basically  the  same  as 
finite  element  shape  functions  (even  for  elements  cut  by  the  crack  surface).  The  EP  shape 
functions  are  much  more  computationally  efficient  than,  e.g.,  moving  least  square  functions 
which  are  used  in  several  meshless  methods.  Moreover,  all  the  computations  can  be  carried 
out  at  the  element  level  as  in  standard  finite  element  codes.  And,  importantly,  the  numer¬ 
ical  integration  of  the  stiffness  and  mass  matrices  can  be  done  as  efficiently  as  in  standard 
finite  element  computations  since  the  intersections  of  the  global  EP  shape  functions  coincide 
with  the  integration  domains.  This  is  in  clear  contrast  with  meshless  methods  based  on  mov¬ 
ing  least  square  functions  where  the  integration  of  the  stiffness  and  mass  matrices  is  very 
computationally  expensive. 

Another  way  of  introducing  a  discontinuity  in  the  EP  shape  functions  is  by  multiplying 
the  partition  of  unity  functions  (pa  by  basis  functions  Lia  that  are  discontinuous.  Prom  the 
definition  of  the  EP  shape  functions  given  in  (3)  it  is  observed  that  if  the  functions  Lia  are 
discontinuous,  the  resulting  EP  shape  functions  are  also  discontinuous.  The  construction  of 
discontinuous  functions  Lia  tfiat  can  approximate  well  the  solution  in  the  neighborhood  of 
a  3-dimensional  crack  front  is  discussed  in  the  next  section.  The  modeling  of  a  crack  in  the 
EPM  is  therefore  done  with  the  aid  of  the  technique  outlined  in  this  section,  which  creates 
discontinuities  in  the  shape  functions  of  elements  cut  by  a  crack,  and  with  the  aid  of  singular 
functions  near  the  crack  front,  which,  in  addition  of  rendering  a  discontinuous  field  at  the 
crack  front,  also  allows  accurate  approximation  of  the  solution  near  the  crack  front  without 
any  mesh  modifications. 
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3  Customized  Shape  Functions  for  an  Edge  in  3D 


There  are  many  classes  of  problems  for  which  the  structure  of  the  underlying  partial  differential 
equation  can  be  exploited.  Oden  and  Duarte  [37,39]  have  demonstrated  how  knowledge  of 
the  solution  of  the  elasticity  equations  near  a  corner  in  two  dimensional  space  can  be  used 
in  a  partition  of  unity  method  to  efficiently  model  the  singularities  that  occur  in  this  class 
of  problems.  This  type  of  singularity  is  resolved  very  poorly  by  polynomial  functions  such 
as  are  used  in  traditional  finite  element  methods,  unless  a  very  refined  mesh  is  used.  In  this 
section,  the  formulation  proposed  by  Oden  and  Duarte  [37,39]  is  extended  to  the  case  of  edges 
in  three-dimensional  space.  Numerical  examples  are  presented  in  Section  13. 

Consider  an  straight  edge  in  three-dimensional  space  as  depicted  in  Fig.  4.  In  the 
figure,  27r  —  a  is  the  opening  angle.  For  example,  in  the  case  of  a  crack,  a  —  2tt.  Associated 
with  the  edge/crack,  there  is  a  Cartesian  local  coordinate  system  (^,7?,  C)  and  a  cylindrical 
coordinate  system  (r,  9,  CO  with  origins  at  (Ox,  Oy,  Oz). 

The  displacement  field  «(r,  9,  C)  in  the  neighborhood  of  the  edge  (for  points  far  from 
its  vertices)  can  be  written  as  [46, 47] 


u{r,  9,  CO 


r  u^{r,9)  ' 

M 

.+4>. 

{  Urf(r,9) 

=  2: 

Af, 

[  u^{r,9)  J 

i=i 

0 

0 

1 

'  %k{r,9)  '' 

> 

-1-  < 

u^{r,9) 

.  uf-{r,9)  J 

(4) 


where  (r,  9,  CO  are  the  cilyndrical  coordinates  relatives  to  the  system  shown  in  Fig.  4,  u^{r,  9), 
Ur,{r,9)  and  u^{r,9)  are  Cartesian  components  of  u  in  the  C~)  and  C—  directions,  respec¬ 
tively,  u^{r,  9),  u^{r,  9)  and  u^{r,  9)  are  functions  smoother  than  any  term  in  the  sum. 


Figure  4:  Coordinate  systems  associated  with  an  edge  in  3D  space. 
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Assuming  that  the  boundary  is  traction-free  and  neglecting  body  forces,  the  functions 

[46,47] 


4y(r,  ^) 


4?  (^)  ^) 


^  { [«  -  Qf  (Af  +  1)]  sin  xfe  -  Af  sin  (Af  -  2)0} 
^  {[k  +  0<‘’(Af  +  1)]  sin  A<‘>«  +  Af  sin  (Af  -  2)#} 


{[k  +  Qf  (Af  +  1)]  cos  Af  S  +  Af  cos  (Af  -  2)«} 


where  the  eigenvalues  a5^\  A^^  are  found  by  solving 


sin  A^^q;  +  Xj  ^  sin  a  =  0 
sin  xf^a  —  A^^^  sin  a  =  0 


In  the  case  of  a  crack,  where  a  =  27r,  the  eigenvalues  are  aJ-^^  =  xf^  =  Aj 


J  +  1 
2 


i>2 


The  material  constant  k  and  G  are 

k  =  3-4.-  ‘5=2(1^ 

where  E  is  the  Young’s  modulus  and  u  is  the  Poisson’s  ratio. 


The  parameters  Qf^  and  Qf^  are  given  by 

(1)  ^  cos(Af  ^  -  1)q;/2  _ 

^  cos(Aj^^  +  l)(x/2  ^ 


sin(Aj^^  —  l)a!/2 
sin(Aj-^^  -I-  1)q:/2 


sin(Af ^  -  1)q:/2  _  .  (2)  cos(Af^  -  1)q!/2 


sin(Aj-  -t-  1)o!/2 


cos(A^^'  +  1)q:/2 


where 


In  the  case  of  a  crack. 


(s)  _  '^j 


Af  - 


Af  + 1 


n(i)_/~^  y  =  3,5,7,... 

I  -A?)  y  =  l,2,4,6,.. 


s  =  l,2 


-1  y=:  1,2, 4,6,... 

-Af  y=3,5,7,... 
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The  functions  ufj  are  given  by  [46]  (assuming  that  the  boundary  is  traction-free  and 
neglecting  body  forces) 


“Cj  -  ' 


r^smXf^e  j  =  1,3,5,... 


2G 

A3) 


[  ^cosXf^e  j  =  2,4,6, 


where 


Af  = 

■’  a 


Before  using  the  above  functions  to  build  customized  element  partition  shape  functions, 
they  have  first  to  be  transformed  to  the  physical  coordinates  (x,  y,  z)  as  follows: 

Define 


s  =  l,2  3  =  1,.. .,M 


Tr'  :  o 


r  ] 

VFT?  1 

►  =  < 

arctan(|)  ► 

Ic' J 

1 

i  c  J 

The  coordinate  system  (^,  rj,  Q  is  shown  in  Fig.  4. 

Next  define 

ul"^{x,y,z)  =  oT 2^ {x,y,z) 

u^^]{x,y,z)  =  oT^^{x,y,z)  s  =  l,2  j  =  l,...,M 

uf) {x,y,z)  =  uf^  oT2^{x,y,z) 

:  {x,y,z)  ^  {^,ri,0 


f  M 

{  X  -  Ox  ' 

1  y  ~  ^ 

IcJ 

[  ^-0.  J 

(5) 

(6) 


(7) 

(8) 


where  e  xR^,  with  rows  given  by  the  base  vectors  of  the  coordinate  system  {^,t],Q 
written  with  respect  of  the  base  vectors  of  the  coordinate  system  {x,  y,  z),  and  O  =  {Ox,  Oy,  Oz) 
are  the  coordinates  of  an  arbitrary  point  along  the  edge  (cf.  Fig.  4).  In  the  above,  u^^]{x,  y,  z). 
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u^^\x,y,z),  s  =  1,2  and  uf){x,y,z)  are  the  components  of  the  displacement  vectors  in  the 
directions  r]  and  (,  respectively,  written  in  terms  of  the  physical  coordinates  x,  y  and  5;.  To 
get  the  components  of  u  in  the  directions  x,  y  and  2  the  following  transformation  needs  to  be 
applied: 


where  R2  =  (R2 


-1\T 


U 


{x) 


'W 


s  =  1, 2  j  =  1, M 


The  construction  of  customized  EP  shape  functions  using  singular  functions  as  defined 
above,  follows  the  same  approach  as  in  the  case  polynomial  type  shape  functions.  The  singular 
functions  are  multiplied  by  the  partition  of  unity  functions  (pa  associated  with  nodes  near  an 
edge/crack.  In  the  computations  of  Section  13,  the  following  singular  functions  are  used  in 
the  construction  of  customized  EP  shape  functions 


^zl  5 


,,(2) 

^xl  5 


U, 


(2) 
yl  5 


(3) 

Uz2 


which  are  given  by 


(2) 

4?  J 


R2 


,-,(1) 


I  Hi  “a 
J  ,-5(1)  f,(2) 

^  “a  "a 


Cl  Hi 


The  customized  EP  shape  functions  are  then  built  as 

^oc  X  {mS  1 >  “i?  >  “i?  >  4?  ’  “S } 

Here,  a  in  the  index  of  a  finite  element  vertex  node  on  or  near  an  edge/crack  in  3D. 


(9) 


4  A  Posteriori  Error  Estimation  for  The  EPM 

Error  estimation  in  PHLEXcrack^'^  is  performed  using  an  explicit  residual  technique.  This 
method  combines  traditional  residual  error  estimation  techniques  with  methods  based  on  in¬ 
terpolation  theory.  The  resulting  formulation  for  the  error  estimate  can  be  written 


N  r 


W  -  £  2iVC|  E  Cf 


k=l 


R 


L2i^k) 


+  c^ 


'^llL2(rtnafifc) 


where  IHI^;  is  the  norm  associated  with  the  bilinear  form,  N  is  the  number  of  domain  subdi¬ 
visions,  Ft  is  the  portion  of  dO.  with  specified  traction  boundary  conditions  and  R  and  f  are 
the  cloud-weighted  interior  and  boundary  residuals 
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r^  =  r^Yl  [(Ma] 

ijJct 

with  (hp)a  —  ha/ Pa  for  points  inside  the  cloud  uja  and  zero  elsewhere. 

Without  determining  the  constants  Ci,  C2  and  C3  we  use  the  error  indicators  ^  as  the 
basis  of  adaptation  and  for  determining  convergence  rates. 


<t>l  = 


R 


L>2{0.k) 


+  Iki 


2 

L2(rtn5nfc) 


5  Extraction  of  Stress  Intensity  Factors: 

The  Least  Squares  Fit  Method 

Once  the  solution  is  obtained  at  some  time  step,  the  amount  and  direction  of  crack  propagation 
over  the  next  time  increment  can  be  determined.  The  crack  front  is  represented  as  a  series  of 
straight  line  segments  connected  at  vertices.  The  stress  intensity  factors  (SIF’s)  are  calculated 
at  the  vertex  points  along  the  crack  front.  Figure  5  represents  the  flowchart  of  the  fracture 
dynamics  process  as  implemented  in  PHLEXcrack^^. 

The  Least  Squares  Fit  method  is  used  for  the  calculations  of  SIF’s.  In  this  method, 
the  SIF’s  are  obtained  by  minimizing  the  errors  among  the  discretized  stresses  calculated  from 
the  solution  and  the  asymptotic  values.  The  method  has  produced  accurate  results  in  finite 
element  settings.  In  this  project,  it  has  been  extended  to  be  used  with  our  three-dimensional 
dynamic  CEPM  model. 

Define  the  least-squares  functional  as: 

:=  (or'*  -0-)y  M  =  /-///  l  = 

where  K'‘j^  is  the  Ith.  stress  intensity  factor  associated  with  mode  M.  The  inner  product  (  ,  )y 
is  defined  as: 

N  j  m  m  \ 

(w,«)y  =  E  i'E'E^iMDrjVjiXa)  Wa{y), 
a-l  Vj=li=l  / 


and  (see  Figure  6) 

u  =  {ui, U2, •  •  • , Um}, V  =  {vi,V2,...,Vm}  are  any  two  vectors  in Je”*, 
<r'*(a;)  is  the  discretized  stress  vector, 

tr{x)  :■=  is  the  asymptotic  stress  vector, 
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PHLEXcrack 


Figure  5:  Flowchart  of  PHLEXcrack^^. 
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Fm(®)  —  ^m(  r)g^Af{9)  are  the  asymptotic  functions, 
y  are  the  coordinates  of  a  vertex  on  the  crack  front, 

Xa  are  the  coordinates  of  a  sampling  point, 

N  is  the  number  of  sampling  points  in  a  domain  centered  at  y, 

D~j  is  the  inverse  of  an  auxiliary  matrix  D.  In  the  current  work,  D  can  be 

selected  to  be  either  the  material  stiffness  matrix  or  the  identity  matrix,  and 

Wa(y)  is  a  weighting  function  associated  with  sampling  point  Xa  given  by 

Wa(y)  =  M - - — yr~  P  typically  equal  to  3-6. 

||y  —  Xgt  Ijjn 


sample  point 


Figure  6:  Domain  used  for  the  Least  Square  Fit  method. 
iC^s  are  found  by  minimizing  the  least-squares  functional: 


dJ 

dK^M 


=  0 


M  =  I-III  /  =  1,1..., /max 


This  leads  to  the  following  system  of  equations: 

III  Qmax 

E  E  KUK,  F'MOy  =  1=1,...,  U. 

M=I  q=l 
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i 


If  only  the  first  terms  (I  =  1)  of  the  three  modes  are  used,  the  three  corresponding 
are  found  by  solving: 


r  (F}.  *'})>' 

L  (Fi„,F», 


(F},Fj,)y 

(F}j,  F}jr)y 


(F/,  F}jj)y 

K}  1 

(0-‘.Fi)y 

(Fjjr,  F}jj)y 

>  =  i 

Ki  \ 

i  (0’*.F)„)y  J 

The  domain  used  for  the  Least  Squares  Fit  method  is  a  cylinder  centered  at  the  vertex 
with  its  axis  along  the  tangent  of  the  crack  front  at  that  vertex.  The  dimensions  of  the  cylinder 
and  the  number  of  sampling  points  in  the  r,  0,  and  2:  directions  are  input  data.  In  addition, 
one  can  choose  the  type  of  the  D  matrix  and  the  weight  p. 

In  addition  to  the  extraction  method  presented  here  and  implemented  in  the  current 
work,  the  user  can  develop  and  add  his  or  her  own  method.  This  will  be  explained  in  Section  11. 


6  Crack  Evolution  Models 

The  Stress  Intensity  Factors  (SIF’s)  calculated  above  are  used  to  determine  whether  the  crack 
will  advance  or  not,  and  the  amount  and  direction  of  propagation,  if  any.  The  front  is  then 
advanced  to  its  new  position,  and  the  numerical  model  updated  accordingly. 

Crack  propagation  quantities  are  calculated  based  on  some  physical  models.  Crack 
physics  are  not  well  known,  particularly  for  three  dimensional  problems.  Therefore,  fracture 
models  usually  make  extensive  use  of  plane  strain  physics  models.  In  this  work,  two  physical 
models  have  been  used.  Moreover,  the  user  can  implement  his/her  own  model  as  described  in 
Section  11. 


The  Freund  Model  [20,24,44]; 

Direction  of  crack  growth  in  the  plane  normal  to  the  crack  front  is  given  by: 


0  = 


2  tan 


(10) 


for  Kii  ^  0,  and  0  =  0  for  Kn  =  0.  In  the  equation  above,  9  is  measured  with  respect  to 
the  forward  vector  ni.  Vector  nj  is  the  crack  front  normal  vector  along  the  intersection  of  the 
plane  normal  to  crack  front  with  the  plane  tangent  to  it  at  the  vertex  (see  Figure  7). 

It  was  assumed  that  mode  III  does  not  affect  crack  direction;  it  only  affects  crack  speed. 
The  sign  in  equation  (10)  above  is  chosen  to  correspond  to  a  positive  stress  intensity  factor 
for  the  direction  given  by  9  (positive  if  Kjj  <  0  and  negative  if  Ku  >  0).  This  corresponds 
to  the  direction  normal  to  the  maximum  hoop  stress. 
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®j  forward  vector 
02  crack  surface  normal 
t  tangent  vector 


Figure  7:  Local  unit  vectors  at  crack  vertex. 


The  current  energy  release  rate  of  a  stationary  crack  is  calculated  at  every  vertex  as: 


G(0) 


1^2 

-^I^equiv 

E* 


+ 


Khi 

2n  ’ 


where  ^ 

Ki,equiv  =  Kicos^{e/2)  -  -itT// cos(0/2)  siu 


H  is  the  shear  modulus,  and  E*  is  the  effective  Young’s  modulus.  For  plane  strain,  E*  is  given 


by: 


E* 


E 

1  —  i/2‘ 


Crack  will  propagate  at  any  vertex  if: 


(?(0)  >  (11) 

where  Ga-it  is  the  critical  energy  release  rate  given  by: 

and  Kid  is  the  dynamic  fracture  toughness  in  a  pure  mode  I  crack;  it  is  approximated  by  the 
constant  value  Kid  (material  property). 

The  speed  at  which  the  crack  will  propagate  at  a  vertex  (d)  is  then  calculated  as  the 
root  of  the  quadratic  equation: 

^^2  _  (ijr  +  i)e  +  =  0 
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where  T  =  G'(0)/Gcrit  >  1,  ^  =  CRjcum  >  1,  and  c  =  d/cR  <  1.  is  the  limiting  crack 
speed  (<  Cr)  given  as  an  input  quantity,  cr  is  Rayleigh  wave  speed  given  as  a  root  of: 


4^/52  ~  (1  + 

where, 

Pl  =  l-  (q  =  dilatational  wave  speed  =  in  piane  strain), 

^1  =  1-  {cs  =  shear  wave  speed  =  in  plane  strain),  and 

K  =  Kosolov  constant  =  3  -  4z/  (for  plane  strain). 

In  the  above,  E  is  Young’s  modulus,  //  is  the  shear  modulus,  p  is  the  mass  density,  and  v  is 
Poisson’s  ratio. 


Prescribed  Velocity  Model; 

In  this  model,  the  propagation  angle  is  calculated  using  equation  (10)  above,  and 
the  propagation  criteria  is  given  by  equation  (11).  The  crack  speed,  however,  is  calculated 
according  to  a  pre-defined  function.  This  function  is  given  as  an  input  data.  The  function  is 
usually  a  segmental  linear  function  of  time  (refer  to  the  User  Manual). 


7  Dynamic  Flow- Structure  Interaction  Problems 

7.1  Introduction 

Although  the  scope  of  the  project  was  restricted  to  the  analysis  of  the  solid  structure  in  the 
dynamic  regime  (with  the  loads  originated  from  the  fluid-structure  interaction),  in  this  section 
we  will  present  short  overview  of  modern  techniques  used  in  this  field.  The  present  section 
will  only  provide  a  brief  overview  of  the  current  state  of  the  art  and  main  contributions  of  the 
Computational  Mechanics  Company  to  these  topics. 

Computational  Fluid  Dynamics  is  presently  dominated  by  finite  volume  and  finite  el¬ 
ement  techniques,  which  provide  greater  flexibility  of  geometry  and  mesh  design  than  earlier 
finite  difference  methods.  There  exist  a  great  number  of  general  purpose  and  specialized  CFD 
codes,  covering  flow  ranges  from  low-speed  non-Newtonian  to  hypersonic  flow  regimes.  Their 
itemization  is  beyond  the  scope  of  this  proposal,  comprehensive  lists  can  be  found  on  the  World 
Wide  Web  (http : //www . comco . com/f eaworld,  http : //icemcf d . com/cf d/CFD_codes . html) . 

COMCO  has  made  significant  and  unique  contributions  to  the  development  of  finite  ele¬ 
ment  methods  for  CFD  applications  including,  in  particular,  adaptive  computational  methods 
for  CFD  [23,36,40-43,53,54].  In  fact,  COMCO  was  the  developer  of  the  first  hp-adaptive  com¬ 
mercial  CFD  code  P3/CFD  marketed  by  PDA  engineering  in  the  early  1990’s.  More  recently. 
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COMCO  has  developed  a  new,  advanced  hp-adaptive  CFD  code  PHLEXcfd  and  a  specialized 
commercial  extrusion  modeling  software  HyperXtrude,  marketed  by  Altair  Engineering. 

Flow-Structure  Interaction  (FSI)  has  been  the  topic  of  many  research  efforts,  especially 
in  application  to  hypersonics  [49]  and  aeroelasticity  [1,8,55].  Historically,  most  of  the  software 
development  efforts  in  this  area  have  been  focused  on  selected  applications  and,  accordingly, 
utilized  some  specific  flow-structure  coupling  schemes  (one-way,  staggered  one-way,  two-way, 
etc.).  To-date,  virtually  the  only  attempt  at  releasing  a  commercial  general  purpose  mul¬ 
tiphysics  code  is  Spectrum  (tm)  by  Centric  -  due  to  very  complex  problem  setup  and  poor 
computational  efficiency,  it’s  success  has  been  rather  limited.  It  appears  that  developing  spe¬ 
cialized  applications  for  specific  types  of  flow-structure  interaction  is  a  more  successful  and 
efficient  approach  than  general  purpose  applications. 

In  the  area  of  flow  interaction  with  deformable  structures,  COMCO  has  developed  a 
unique  module  for  modeling  closely  coupled  dynamic  flow-structure  interactions.  In  particular, 
this  module  utilizes: 

•  large  deflection  structural  model, 

•  arbitrary  Lagrangian-Eulerian  fluid  solver,  and 

•  a  moving  mesh  algorithm, 

to  model  fully  coupled  dynamic  flow-structure  interactions  representative  of  flutter  and  buffet¬ 
ing.  This  application  has  been  developed  under  sponsorship  from  Air  Force  Office  of  Scientific 
Research  (AFOSR)  [55]. 

COMCO  has  been  developing  computational  software  for  various  government  installa¬ 
tions  and  private  industry  for  over  a  decade.  A  few  of  these  development  efforts  which  are 
relevant  are  as  follows: 

•  -  Title:  Adaptive  Computational  Methods  for  Fluid  Structure  Interaction  in  Internal 

Flows 

-  Description:  This  project  focused  on  the  development  of  an  inviscid  h-adaptive 
finite  element  compressible  flow  solver  for  modeling  the  transonic  and  supersonic 
flow  through  a  jet  engine.  In  particular,  this  effort  used  mesh  sliding  interface 
capabilities  to  move  the  rotor  blade  relative  to  the  stator. 

-  Client:  NASA,  contract  NAS3-25196. 

-  Contact:  Jose  Sans 

-  Date  of  Completion:  May  1989. 

•  -  Title;  Analysis  of  Flow-,  Thermal-  and  Structural-Interaction  of  Hypersonic  Struc¬ 

tures  Subjected  to  Severe  Aerodynamic  Heating 
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—  Description:  The  project  was  focused  on  development  of  finite  element  techniques 
for  coupling  transient  hypersonic  flow  simulation  and  heating  rate  prediction  with 
thermal-structural  response  of  hypersonic  structures. 

—  Client:  AFOSR 

—  Contact:  Dr.  Spencer  Wu,  (phone:  703-696-8523) 

—  Date  of  Completion:  May  1990. 

•  -  Title:  STORESIM:  An  Integrated  System  for  Multi-Body  CFD  Simulations  Using 

Unstructured,  Adaptive  Grids. 

-  Description:  The  primary  objective  of  this  project  was  to  develop  a  state-of-the 
art  integrated  software  package  to  accurately  and  efficiently  model  the  flow  around 
the  missiles,  predict  the  trajectories  by  incorporating  the  fluid-structure  interaction 
during  missile  deployment. 

—  Client:  Wright-Patterson  Air  Force  Base 

-  Contact:  Frank  C.  Witzeman  (phone:  937-255-3778,  e-mail:  witzeman@ind4.fim.wpafb.af.mir 

—  Date  of  Completion:  December  1996. 

•  -  Title:  Dynamic  Instabilities  in  Flow-Structure  Interactions,  Flutter  and  Buffeting 

-  Description:  The  project  was  focused  on  a  detailed  study  of  dynamic  instabilities 
in  fluid-structure  interaction  and  development  of  new  computational  techniques 
for  the  numerical  simulation  of  a  class  of  physical  instabilities  in  which  a  viscous 
compressible  fluid  interacts  with  a  flexible  elastic  structure  undergoing  large  defor¬ 
mations. 

-  Client:  AFOSR 

-  Contact:  Mjr.  Brian  Sanders  (phone:  703-696-7259,  e-mail:  brian.sanders@afosr.af.mil) 

—  Date  of  Completion:  December  1997. 

The  basics  of  linear  flow-structure  interactions  have  been  covered  in  many  monographs 
(see  for  example  [7, 12])  and  numerous  scientific  papers.  Presently,  most  of  the  methods  in 
applied  aeroleasticity  rely  strongly  on  linearization  of  the  equations  of  motion  and  simplified 
computational  methods,  such  as  analytical  solutions,  vortex  lattice  and  doublet  lattice  models, 
or  finite  volume  discretizations. 

In  the  design  of  modern  aerospace  systems,  the  necessity  of  considering  more  complex, 
nonlinear  physical  phenomena  in  fluid  structure  interactions,  such  as  those  encountered  in 
buffeting,  fluid  thermal  structural  interactions,  and  viscous  turbulent  boundary  layer  effects 
has  pushed  the  analysis  problem  far  beyond  the  limits  of  the  classical  theory  and  into  a  realm 
of  phenomena  for  which  a  list  of  open  physical,  mathematical  and  computational  issues  exist. 

These  problems  require  a  fully  nonlinear  formulation  of  flow-structure  interaction  problems,  a 
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more  general  notion  of  stability,  applicable  to  nonlinear  equations  characterizing  viscous  flow 
interacting  with  a  complex  structure  undergoing  large  deformations,  and,  possibly,  large  tran¬ 
sient  thermal  gradients.  The  linearized  stability  elements  of  the  classical  theory  are  inadequate 
and  new  approaches  to  these  classes  of  problems  are  needed. 

Recently,  several  efforts  have  been  reported  towards  expanding  the  flow-structure  inter¬ 
action  solution  into  the  realm  of  nonlinear  flow  analysis  coupled  in  some  way  with  structural 
deformations.  Bendiksen  [6]  coupled  a  finite  difference  Euler  CFD  scheme  with  a  simplified 
finite  element  structural  representation  to  solve  transient  flow-structure  interaction  problems. 
Alonso  et.al  [1],  Blackburn  et.al.  [8]  and  Morton  et.al.  [34]  studied  experimentally  and  numer¬ 
ically  the  unstable  oscillation  of  a  rigid  cylinder  induced  by  interaction  with  a  uniform  flow. 
Flutter  of  a  nonlinear  buckled  plate  was  studied  by  Dowell  [11]  and  Virgin  and  Dowell  [56]. 
Various  nonlinear  aspects  of  aeroelasticity  were  discussed  by  Dugundji  [18]. 

7.2  A  General  Problem  Statement 

This  section  presents  a  detailed  formulation  of  dynamic  flow-structure  interaction  problems 
and  computational  techniques  used  to  solve  them.  This  formulation  was  developed  during  the 
“Flutter  and  Buffeting”  project  mentioned  above: 

•  it  is  applicable  to  fully  coupled  flow-structure  interaction  problems,  with  complete  preser¬ 
vation  of  dynamic  coupling  terms  on  the  interface, 

•  it  can  correctly  represent  large  structural  deflections  and  their  interaction  with  the  flow- 
field, 

•  it  captures  the  nonlinearities  that  contribute  to  instabilities  and  define  the  final  motion 
of  the  system  (such  as  limit  cycle  oscillations),  and 

•  it  lends  itself  to  finite  element  discretization  or  other  numerical  approximation  methods. 

The  formulation  combines  an  Arbitrary  Lagrangian-Eulerian  (ALE)  approach  for  the 
fluid  motion  with  an  Updated  Lagrangian  Reference  (ULR)  description  for  the  solid.  The 
general  idea  and  details  of  this  formulation  are  presented  in  the  reminder  of  this  section. 

A  general  illustration  of  the  problem  under  consideration  is  shown  in  Fig.  8.  Consider  a 
solid  body  with  a  stress-free  configuration  Qq.  At  the  time  t,  the  solid  occupies  a  configuration 

In  this  configuration  the  solid  is  interacting  with  fluid,  with  a  finite  element  or  other  mesh 
discretizing  the  flow  domain.  Because  of  the  flow-structure  interaction  or  other  forces,  the  solid 
is  moving  towards  the  new  current  configuration  Qt+dt-  Accordingly,  the  mesh  discretizing  the 
fluid  domain  is  also  moving  so  that  the  boundary  of  the  flow  domain  continues  to  match  the 
new  position  of  the  solid.  The  mesh  position  is  updated  at  every  time  step. 
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Figure  8:  A  general  statement  of  the  flow-structure  interaction  problem. 

The  details  of  the  corresponding  equations  for  the  solid,  fluid  and  the  coupled  problem 
are  presented  below. 


7.3  Structural  Deformation  Problem 

This  section  presents  the  formulation  of  the  structural  deformation  problem.  The  main  objec¬ 
tive  is  to  develop  a  formulation  that  is  applicable  to  large  deflection  structural  problems  and 
is  convienient  for  close  coupling  with  flow  equations  and  moving  computational  meshes. 

The  main  assumptions  considered  for  the  structural  problem  are  as  follows: 


•  large  deformations  and  small  strains,  and 

•  elastic  material  behavior. 

The  approach  to  the  description  of  this  problem  is  illustrated  in  Fig.  9.  All  the 
positions  and  displacements  are  represented  in  an  inertial  coordinate  system  Xj,  i  =  1,2,3. 
At  some  initial  time  r  =  to  the  body  occupies  an  initial,  stress-free  configuration  For  a 
particle  P,  the  position  vector  in  this  configuration  is  defined  by  a  vector  Xq.  At  a  certain 
time  r  =  i  the  body  is  in  a  configuration  Slj,  and  a  position  vector  for  a  particle  P  is  defined 
as  X.  Typically,  in  the  present  solution  process,  this  is  the  configuration  reached  in  a  transient 
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Updated  Lagrangian 


Figure  9:  Transient  structural  deformation  at  large  deformations. 

solution  process.  At  the  most  current  time  t  =  t  +  dt  the  body  is  in  a  configuration 
with  the  position  vector  for  a  particle  P  defined  as  r.  Typically,  this  is  the  configuration  being 
calculated  in  a  transient  solution  process. 

Generally,  the  deformation  can  be  viewed  from  an  arbitrary  reference  configuration  fi, 
which  need  not  coincide  with  any  of  the  positions  of  the  body  during  the  motion.  In  the  present 
flow-structure  interaction  problem,  the  reference  configuration  is  chosen  to  be  wherever  the 
computational  mesh  is  located.  Typically,  this  will  be  the  latest  converged  configuration  O*, 
as  shown  in  Fig.  9.  Importantly,  all  the  differentiation  and  integration  discussed  subsequently 
will  be  performed  on  this  reference  configuration. 

It  should  be  noted  that  the  reference  configuration  is  a  Lagrangian  reference,  serving 
as  a  convenient  viewpoint  on  the  total  deformation  from  fio  to  Q.t+dt-  As  the  solution  pro¬ 
gresses  and  the  mesh  moves,  different  configurations  will  be  chosen  as  this  reference  frame. 
Thus,  the  present  formulation  was  named  an  Updated  Lagrangian  Reference  (ULR)  formula¬ 
tion.  Importantly,  it  differs  from  the  Updated  Lagrangian  formulation  in  this  that  the  total 
strains  and  stresses  do  not  get  accumulated  during  the  computational  process,  but  are  always 
recalculated  from  the  total  deformation.  Thus,  accumulation  of  errors  is  avoided. 

Within  the  above  framework,  the  solution  of  a  transient  deformation  problem  can  be 
stated  as:  Knowing  the  present  updated  Lagrangian  reference  configuration  and  the  loads 
acting  on  the  body,  calculate  the  new  current  configuration  Q^t+dt  after  a  time  interval  dt. 
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Equations  of  Motion.  Starting  at  this  point,  in  order  to  differenciate  between  the  solid  and 
fluid  domain,  the  domain  corresponding  to  the  solid  body  will  be  indicated  by  a  superscript 
b.  The  equations  of  motion  for  the  structure,  expressed  in  the  reference  configuration  f2^,  are 
of  the  form: 


Pb  a{T)  =  diva:  T(r)  +  ps  b('r)  in  OP  x  [0,  T]  (12) 

where  ps  is  the  density  of  the  body  in  the  reference  configuration,  r  is  the  time,  a  is  an 
acceleration  vector,  T  is  a  first  Piola-Kirchhoff  stress  tensor,  6  is  a  body  force  vector  per  unit 
mass  and  diVa,  is  a  divergence  operator  with  respect  to  the  reference  configuration. 

The  equations  of  motion  can  also  be  expressed  in  terms  of  the  second,  symmetric 
Piola-Kirchhoff  stress  tensor  as: 

Pb  a  =  diva,  {FS)  ■+•  pn  b  (13) 

where  the  time  symbol  was  dropped  for  clarity,  S  denotes  the  second  Piola  Kirchhoff  stress 
tensor  and  jP  is  a  deformation  gradient  with  respect  to  the  reference  configuration  Q®: 

F  =  grada;r  (14) 

Constitutive  Equations.  For  the  materials  considered  in  this  project,  the  constitutive 
equations  are  defined  by  a  linear  relationship  between  the  Cauchy  stress  and  strain  tensors. 
It  can  be  shown,  using  standard  arguments  of  Continuum  Mechanics  [22,49],  that  for  large 
deflections  /  small  strains  problems  this  is  equivalent  to  a  linear  relationship  between  the 
second  Piola-Kirchhoff  stress  tensor  and  a  Green-Lagrange  strain  tensor: 


S  =  CE  (15) 

where  C  is  a  fourth  order  elasticity  tensor,  satisfying  the  standard  conditions  of  boundedness, 
symmetry  and  ellipticity.  For  isotropic  materials,  this  can  be  expressed  as: 

S  =  \tr{E)I  +  2pE  (16) 

where  A  and  p  are  Lame  constants  and  I  is  an  identity  tensor. 


Kinematics  and  Strain  Measures.  In  classical  continuum  mechanics,  the  Green  strain 
tensor  is  defined  as: 
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(17) 


E  =  i  (F^F  -  I) 

where  F  is  a  deformation  gradient  from  the  stress-free  configuration  to  the  deformed  config¬ 
uration.  This  formula  is  valid  under  the  usual  assumption  that  the  reference  configuration 
coincides  with  the  initial,  stress-free  configuration. 

However,  for  the  ULR  description  introduced  above  this  is  not  the  case,  and  a  more  geral 
formula  must  be  used.  Considering  that  the  Green  strain  tensor  measures  the  change  of  metrics 
between  the  initial  and  final  configuration  as  “seen”  from  the  reference  configuration  [52],  a 
generalized  formula  for  the  strain  tensor  in  the  ULR  description  is  defined  as: 

E  =  ^  {F'^F  -  Fo'^Fo)  (18) 

where  Fq  is  the  deformation  gradient  from  the  reference  configuration  to  the  initial,  stress-free 
configuration: 


Fq  =  grad^  xq  (19) 

It  can  be  easily  verified  that  the  formula  (18)  reduces  to  the  Green  strain  tensor  or 
Almansi-Hamel  strain  tensor  if  the  reference  configuration  coincides  with  the  initial  or  the 
deformed  configurations,  respectively. 


Boundary  Conditions.  For  the  above  problem,  the  typical  boundary  conditions  include: 
the  prescribed  motion  on  a  part  of  the  boundary,  presecribed  traction  history,  or  their  combi¬ 
nations.  The  prescribed  motion  is  defined  as: 

r(x,  r)  =  r{x,T)  on  x  [0,T]  (20) 

where  Fu  is  the  portion  of  the  boundary  with  a  prescribed  displacement  history  and  f  is  the 
prescribed  motion  vector. 

The  boundary  tractions  t{x,  r)  on  the  solid  surface  are  expressed  in  terms  of  the  first 
Piola-Kirchhoff  stress  tensor  as: 


t  =  Tn  (21) 

where  n  is  a  vector  normal  to  the  boundary  in  the  reference  configuration.  Since  this  is  the 
configuration  discretized  by  a  finite  element  mesh,  the  calculation  of  boundary  normals  is 
straightforward.  For  a  portion  of  the  boundary  Ft,  a  traction  history  is  prescribed  as: 
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t(x,T)  =  i{x,T) 


on  Ft  X  [0,T] 


(22) 


where  i  is  the  prescribed  traction  vector.  On  the  flow-structure  interface  F^,  the  tractions 
on  the  solid  surface  are  a  result  of  fluid  pressure  and  viscous  shear.  The  details  of  these 
interactions  are  presented  later  in  this  document. 


Weak  Formulation.  Using  a  standard  approach,  namely  multiplying  (13)  by  a  test  function 
and  integrating  over  the  reference  domain  the  following  weak  formulation  of  the  equations 
of  motion  is  obtained: 

Find  a  mapping  r(r)  from  [0,T]  to  V  such  that: 

f  pbw  a  dQ+  f  (gradj-u;)^  :  (FS)  dQ  =  (23) 

JnB  JciB 

f  pbW  •  b  dO,  +  f  w  •  i  dr  +  f  w  •  {FS)n^  dT 
Job  Jvt  Jrt 

where  V  =  {w  e  (w)  =  0  on  F„}  is  the  space  of  admissible  displacements  for 

the  structural  domain,  with  N  denoting  the  dimensionality  of  the  space  (iV  =  2  or  3).  In 
the  above  formulation,  the  Dirichlet  boundary  conditions  on  F„  are  a  part  of  the  definition 
of  the  space  U  -  in  practice  they  are  enforced  through  a  penalty  method  or  via  direct  matrix 
elimination.  At  this  point,  we  will  leave  the  boundary  terms  on  the  flow-structure  interface  F^ 
in  the  above  generic  form  -  the  handling  of  flow-structure  interactions  will  be  discussed  later. 


7.4  Flow  Equations;  an  ALE  Description 

This  section  presents  flow  equations  in  Arbitrary  Lagrangian-Eulerian  description,  which  uti¬ 
lizes  a  moving  reference  frame  (grid)  to  handle  changes  in  the  flow  domain  due  to  structural 
motion. 


Governing  Equations.  The  transient  incompressible  flow  in  an  ALE  description  is  gov¬ 
erned  by  the  momentum  equations: 


PfU  +  Pf{{u-  u^)  •  grad)u  -diver  =  f  in  OF  x  [0,  T]  (24) 

augmented  by  the  incompressibility  condition: 

div  u  =  0  (25) 
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Here  pp  is  the  density  of  the  fluid,  u  is  the  flow  velocity,  is  a  prescribed  grid  velocity,  cr 
is  the  stress  tensor,  /  represents  body  forces  acting  on  the  fluid,  and  denotes  the  flow 
domain. 

The  grid  velocity  in  the  above  equation  is  a  prescribed  field  u^{x,  r),  which  is  generally 
quite  arbitrary.  In  the  approach  adopted  in  this  work,  the  grid  velocity  on  the  fluid-solid 
interface  is  equal  to  the  velocity  of  the  solid  body,  so  that  the  grid  follows  the  motion  of  the 
body. 


The  stress  tensor  in  the  fluid  is  expressed  as: 


or  =  —pi  +  z/(grad«  -I-  (gradw)^)  (26) 

where  p  is  the  hydrostatic  pressure  and  v  is  fluid  viscosity. 

In  numerical  simulations  it  is  often  convenient  to  relax  the  incompressibility  constraint 
and  introduce  a  penalty  approach  to  enforce  its  approximate  satisfaction.  In  such  cases,  the 
pressure  is  expressed  as: 


p  =  —A  div  u 


(27) 


where  A  is  a  (large)  penalty  number. 

The  typical  boundary  conditions  for  the  flow  domain  include  prescribed  velocity  on  the 

inflow: 


u{x,t)  =  u{x,t)  on  Tj  x  [0,T]  (28) 

prescribed  traction  (usually  dominated  by  pressure): 

t{x,  r)  =  i{x,  t)  onTp  X  [0,  T]  (29) 

or  solid  wall  no-slip  boundary  condition: 

u  —  =  0  on  X  [0,  T].  (30) 


Weak  Formulation.  Quite  similarly  as  for  the  solid  domain,  the  weak  formulation  for  the 
flow  equations  is  obtained  by  multiplying  the  governing  equations  by  a  test  function  and 
integrating  over  the  domain  of  the  fluid.  After  integrating  relevant  terms  by  parts,  the  final 
weak  statement  becomes: 
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Find  a  mapping  {u,p)  from  [0,T]  toV  xQ  such  that: 


{ppiijw)  +  c{u,u,w)  +  a{u,w)  +  b{w,p)  —  {f,w)  +  b.i.  (31) 

b{u,q)  =  0 

for  every  {w,  q)  inV  x  Q.  The  consecutive  terms  are  defined  as: 


{PfU,w)  = 

/  pfU-w  dQj 

c{u,  u,  w)  = 

/  pfu-  grad u-wdQ 

Jaf 

a(u,  w)  = 

[  i/(grad  u  +  (grad  u)^)  :  grad  w  dVt 

b{w,p)  = 

—  /  div  w  p  dQ 

Jn^ 

and  the  boundary  contributions  denoted  as  b.i.  are  expressed  as: 


J  [z/(gradu  +  (gradu)^)n^  -pn^]  • 

over  the  entire  boundary  of  the  flow  domain  (including  the  solid-fluid  interface  F  ).  In  the 
above,  V  and  Q  are  the  spaces  of  admissible  flow  velocities  and  pressures:  V  =  {w  e 
(B^(n))^,  (w)  =  0  on  Ti}  and  Q  =  {q  e  respectively.  The  traction  boundary 

conditions  are  applied  by  prescribing  the  value  of  the  traction  in  the  boundary  integral: 


/  i •  w  dr 

Jr„ 

and  the  Dirichlet  boundary  conditions  on  the  inflow  and  solid  wall  are  usually  enforced  via 
penalty  or  direct  matrix  operations. 


7.5  Flow- Structure  Interaction  Equations 

This  section  addresses  the  formulation  of  the  dynamic  flow-structure  interaction  problem  at 
large  structural  deflections.  The  main  objectives  of  this  formulation  are: 

•  to  facilitate  a  fully  coupled  solution  of  the  structural  and  the  flow  problem  in  a  single 
time  marching  procedure,  and 

•  to  handle  the  dynamic  interaction  between  the  solid  and  the  fluid  with  minimum  com¬ 
putational  cost  or  complexity. 
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a)  Solid  Body  Domain 


b)  Fluid  Domain 


Figure  10:  Fluid-structure  interaction:  complementary  subdomains  and  interface  conditions. 

The  structural  deformation  and  fluid  equations  formulated  so  far  are  essentially  dis¬ 
joint,  in  the  sense  that  they  describe  the  behavior  of  the  fluid  and  the  solid  separately,  without 
any  implication  of  interaction  between  these  media.  However,  several  characteristics  of  these 
formulations  facilitate  their  application  to  the  solution  of  a  fully  coupled  flow-structure  inter¬ 
action  problem,  namely: 

1.  By  selecting  the  reference  configuration  for  the  ULR  structural  description  to  coincide 
with  the  most  recent  configuration  attained  during  the  motion,  the  weak  formulations 
for  the  structural  and  flow  problem  span  complementary  subdomains  (see  Fig.  10): 

=  Q®  U 

Therefore,  a  single  mesh  covering  both  the  fluid  and  solid  computational  domains  can 
be  used  in  numerical  simulations.  Note  that  this  would  not  be  the  case  if  the  total 
Lagrangian  description  were  used  for  the  solid  -  the  fixed  solid  computational  domain 
would  be  “lagging”  behind  the  changing  flow  domain. 

2.  Owing  to  the  ALE  formulation  for  the  fluid  domain,  the  grid  motion  can  be  used  to 
continuously  match  the  deformation  of  the  solid. 

3.  Both  the  structural  and  flow  equations  represent  momentum  balance  in  true  physical 
units.  Therefore,  the  weak  forms  are  compatible  in  the  sense  that  the  terms  of  each 
weak  form  represent  virtual  power  and  can  be  added  together. 

In  the  context  of  interaction  problems,  of  special  interest  are  the  integrals  on  the  flow- 
structure  interface  F^.  It  can  be  easily  observed  that  the  corresponding  boundary  terms  for 
the  structural  problem  represent  surface  tractions  on  the  solid  body  (see  Fig.  10): 

w  ■  {FS)n^  dT  =  dr  (32) 
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and  the  similar  also  holds  for  the  flow  equations: 

J  w  •  [i/(grad  u  +  (grad  w)^)n^  —  pn^]  dT  —  w  cfT.  (33) 

The  interface  conditions  for  the  flow-structure  interaction  problem  require  that  the 
motion  of  the  grid  matches  the  motion  of  the  solid  body: 

=  r  on  (34) 

and  that  the  fluxes  (tractions)  cancel  out  across  the  interface: 

on  (35) 

In  the  context  of  the  weak  statement,  the  total  integral  on  the  flow-structure  interface 

consists  of  both  fluid  and  solid  contributions.  In  the  light  of  equations  (32)  and  (33),  it 
amounts  to:  . 

+  (36) 

which,  because  of  the  interface  condition  (35),  is  equal  to  zero.  This  means  that  in  a  solution 
scheme  that  utilizes  the  formulation  presented  above  to  simultaneously  solve  the  union  of 
a  solid  and  a  fluid  domain,  there  is  no  need  to  evaluate  additional  interface  integrals  -  the 
interface  terms  cancel  out.  This  is  subject  to  additional  conditions,  such  as  correctly  calculated 
grid  motion  and  using  time  integration  schemes  which  do  not  introduce  additional  boundary 
terms  as  a  by-product  of  the  algorithm. 

Also,  it  should  be  noted  that  the  primary  variables  in  the  structural  and  fluid  weak 
statements  are  different  (positions  vs.  velocities),  which  complicates  the  direct  formulation  of 
a  single  problem  encompassing  both  the  fluid  and  the  solid  subdomains.  This  situation  will 
be  rectified  by  selecting  appropriate  time  discretization  schemes,  with  velocity  becoming  the 
primary  variable  for  the  entire  domain.  The  details  are  discussed  in  the  next  section. 


8  Solution  of  the  Flow-Structure  Interaction  Problem 


In  order  to  solve  the  flow-structure  interaction  problem  without  introducing  time  lag  error 
typical  of  staggered  solution  sequences,  the  approach  developed  in  this  project  includes  the 
entire  domain  in  the  computational  scheme  and  handles  the  solid  and  the  flow  equations 
simultaneously.  To  facilitate  this,  appropriate  versions  of  time  marching  algorithms  have  been 
introduced  in  which  velocities  are  the  primary  unknowns  for  both  the  solid  and  fluid  domains. 
In  particular,  in  this  work  we  used: 
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•  for  the  solid  domain,  a  one-step  Newmark  algorithm  in  velocity  formulation,  and 

•  for  the  flow  domain,  a  one-step  ^-method. 

Based  on  the  comments  from  the  previous  section,  we  expect  the  interface  integrals 
to  cancel  out,  to  simplify  the  computational  procedure.  Additionally,  a  mesh  movement 
algorithm  will  be  used  to  continually  update  the  computational  domain.  The  basic  steps  and 
formulas  of  the  numerical  solution  are  discussed  further  in  this  section. 


8.1  Time  Integration  for  the  Structural  Domain 

To  advance  the  solution  in  time  on  the  structural  subdomain,  a  Newmark  family  of  algorithms 
was  used.  The  basic  formulas  of  the  Newmark  family  express  positions  and  velocities  at  the 
time  step  n  -I- 1  as: 

*"+1  =  x”  +  dtu” -H -^[(1  -  2^)a"  + 

-I-  (1  —  a”  +  7  dt 

where  the  index  n  indicates  the  previous,  converged  time  step.  ^  and  7  are  Newmark  parame¬ 
ters,  which  affect  dissipation  and  dispersion  properties  of  the  algorithm.  In  order  to  formulate 
a  velocity  version  of  this  algorithm,  the  above  formulas  were  recast  to  represent  new  positions 
and  accelerations  in  terms  of  velocities: 

j  2  7 

0”+^  =  -  v”)  - 

jdt  7 

Substitution  of  these  formulas  to  the  weak  statement  for  the  solid  leads  to  the  following 
formulation: 


(grad^iy)^  :  dQ  =  (37) 

f  PbW-(-^v”'  +  - — -d^)dCl+f  pbW  ■  dQ+  [  w-i^^^dT 
7dt  7  Ins  Vr* 

This  is  an  equation  for  calculating  the  velocities  at  a  new  time  step  Due  to  the 
implicit  dependence  of  the  nonlinear  terms  upon  the  solution  at  the  new  time  step, 

this  is  actually  a  nonlinear  equation.  It  is  solved  using  a  Newton  iterative  scheme  based  on 
a  linearization  of  equation  (38)  with  respect  to  velocities  The  details  of  this  procedure 
follow  the  standard  approach  for  nonlinear  problems  and  will  not  be  presented  here. 
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8.2  Time  Integration  for  the  Flow  Domain 

To  advance  the  solution  for  the  flow  domain,  a  one-step  implicit  0-method  was  used.  It  is 
based  on  an  approximation  of  the  rates  of  change  of  the  solution  as. 

where  9  e  (0, 1].  It  is  a  first  order  algorithm  except  for  0  =  i,  at  which  it  has  a  second  order 
time  accuracy.  Substitution  of  the  above  formula  into  the  weak  statement  for  the  fluid  leads 

to  the  following  equation: 


w)  +  +  a(«"+\  to)  + 

'  0  dt 

b(u^^\q) 


0 


(39) 


where  the  consecutive  terms  were  explained  in  Section  7.4. 

The  above  equation  calculates  velocities  at  a  new  time  step  Due  to  the 

presence  of  nonlinear  terms  this  is  actually  a  nonlinear  equation.  Similarly 

as  for  the  solid  subdomain,  it  is  is  solved  using  a  Newton  iterative  scheme  based  on  linearization 
with  respect  to  velocities 


8.3  Solution  of  Combined  Flow-Structural  Equations 


A  significant  characteristic  of  the  time-discretized  structural  and  flow  equations  (38)  and  (40) 
is  that  they  calculate  velocities  at  the  new  time  step  T+i.  Thus,  the  coresponding  Newton 
equations  can  be  discretized  and  solved  for  the  entire  domain  simultaneously.  In  particular, 
discretization  by  the  finite  element  shape  functions  leads  to  the  following  incremental  matrix 
equations  to  be  solved  in  each  nonlinear  iteration: 


Ai'^^  dUi'^^  =  in  Q,  (40) 

^  ijn+l  +  ^jjn+l 

where  A  is  the  left-hand  side  matrix,  U  is  the  discrete  solution  vector  and  R  is  a  nonlinear 
residual  vector  calculated  according  to  the  Newton  procedure.  The  discrete  solution  vector 
U  represents  velocities  of  both  the  fluid  and  solid  subdomains.  This  simultaneous  solution 
preserves  the  relavant  two-way  coupling  between  the  fluid  and  the  solid  domain. 
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8.4  Moving  Mesh  Algorithm 


Upon  examination  of  equation  (40)  it  can  be  seen  that  the  grid  or  mesh  velocity  appears 
explicitly  in  this  equation.  During  the  solution  of  iterative  equation  (40),  the  grid  velocity  is 
treated  as  known.  However,  the  distribution  of  mesh  velocity  and  position  need  to  be  updated 
after  each  time  step. 

Importantly,  for  the  sake  of  the  solution  of  the  nonlinear  equation  (40),  the  distribution 
of  the  mesh  velocity  is  quite  arbitrary.  The  only  stringent  requirement  comes  from  the  interface 
condition  (34),  which  indicates  that  the  grid  velocity  on  the  fluid-solid  boundary  is  equal  to 
the  velocity  of  the  solid. 

In  the  solution  process,  the  mesh  motion  is  calculated  after  calculating  the  velocities 
yn+i^  The  typical  mesh  motion  problem  and  boundary  conditions  are  characterized  by  a 
Laplace  equation  for  each  of  the  components  of  mesh  velocity: 

fJ'kni  =  0  inQ^  (41) 

where  i,k  =  l,iV  and  uf  are  the  x,y,z  components  of  the  grid  velocities.  The  boundary 
conditions  are: 

ul  =  nf  on 
=  0  on 
ul,n  =  0  on  r® 

Here  is  the  flow-structure  interface,  is  a  fixed  solid  wall  (if  present)  and  indicates  fluid 
boundary  with  other  boundary  conditions,  such  as  prescribed  inflow,  outflow  or  symmetry. 
The  comma  indicates  a  derivative  and  n  denotes  a  direction  normal  to  the  boundary. 

The  above  approach  allows  for  modeling  of  multiple  bodies  with  arbitrary  surface/interface 
movements  as  one  would  expect  in  fluid-structure  interaction  problems. 


8.5  Overall  Solution  Algorithm 

The  dynamic  flow-structure  interaction  problem  is  solved  via  a  time  stepping  algorithm,  with 
a  relevant  update  to  the  mesh  to  accomodate  structural  motion.  The  overall  solution  sequence 
is  as  follows: 

1.  Define  a  finite  element  mesh  covering  both  the  structural  and  flow  domains  in  some 
initial  configuration.  Identify  the  material  properties  and  boundary  conditions. 

2.  Define  an  initial  flow  solution  by  solving  the  steady-state  flow  equations  or  importing  a 
flow  solution  from  an  independent  flow  solver. 
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3.  Perform  time  stepping  from  time  t  =  0  to  the  final  time  t  =  T  using  time  increments  dt. 

For  each  time  step: 

(a)  calcuate  the  velocities  for  the  entire  domain  by  solving  the  flow-structure  interaction 
problem  (40), 

(b)  based  on  the  calculated  velocities,  determine  new  structural  positions  r, 

(c)  update  the  structural  reference  configuration  Q  to  correspond  to  the  new  structural 
positions,  and 

(d)  calculate  the  mesh  velocity  w®  from  equation  (41)  and  update  the  mesh  location. 

The  above  sequence  can  be  augmented  by  an  additional  perturbation  of  the  steady-state 
equilibrium  configuration  (to  trigger  instabilities) ,  saving  the  solution  to  a  file  or  performing 
postprocessing.  Moreover,  at  selected  time  intervals  the  mesh  adaptation  can  be  performed. 

9  Underwater  Shock  Analysis 

9.1  Introduction 

The  development  of  numerical  techniques  to  assess  damage  to  naval  structures  due  to  under¬ 
water  explosions  is  a  continuing  effort  to  ensure  that  our  naval  forces  are  placed  at  less  risk  in 
combat.  Underwater  shock  damage  typically  combines  the  effects  of  the  incident  shock  wave 
that  propagates  at  high  velocity  outward  from  the  explosive  source  and  envelops  the  structure 
with  later  low  frequency  pulsations  of  the  explosion  gas  bubble.  The  shock  wave  induces  high 
frequency  response  in  the  structure  with  possible  localized  damage  while  the  bubble  oscilla¬ 
tion  excites  low  frequency  beam  modes  of  the  entire  vessel  that  can  lead  to  bend-buckling 
phenomena.  Another  source  of  damage  is  due  to  the  possible  collapse  of  the  gas  bubble  on 
the  structure  resulting  in  fluid  impact.  This  latter  effect  can  be  the  most  severe  result  of  an 
underwater  explosion  if  it  occurs  very  close  to  the  structure. 

Naval  vessels  are  typically  constructed  using  discretely  stiffened  plate  and  shell  modules. 
As  a  result,  the  effect  of  an  underwater  shock  is  to  induce  deformations  in  the  structure  that, 
if  they  are  large  enough,  will  produce  compressive  loads  in  the  stiffeners  causing  the  torsional- 
flexural  instability  called  ’’tripping.”  One  result  of  this  process  is  the  tendency  of  the  stiffener 
to  tear  away  from  the  hull  plating  that  may  or  may  not  open  the  hull  to  flooding.  In  such 
a  scenario  the  study  of  the  combined  effects  of  underwater  shock  and  crack  propagation  is  of 
vital  interest  to  the  naval  community.  The  interfacing  of  two  computer  codes  that  address 
these  phenomena  separately  is  therefore  appropriate  to  study  the  interaction  of  these  combined 
effects.  These  codes  are  the  Underwater  Shock  Analysis  (USA)  code  developed  and  supported 
by  Unique  Software  Applications  of  Colorado  Springs,  Colorado,  and  the  PHLEXcrack  code 


developed  and  supported  by  the  Computational  Mechanics  Company  of  Austin,  Texas.  For 
more  details  on  the  theory  and  application  of  USA  code,  one  can  refer  to  [10]. 

9.2  Underwater  Shock  Analysis  Code 

The  USA  code  determines  the  transient  response  of  a  totally  or  partially  submerged  structure 
due  to  an  incident  acoustic  shock  wave  and  subsequent  explosion  gas  bubble  pulsation.  USA 
is  built  around  the  family  of  Doubly  Asymptotic  Approximations  (DAA)  that  treat  the  fluid 
structure  interaction  by  a  boundary  element  approach  that  avoids  the  incorporation  of  a  large 
volume  of  fluid  around  the  structure  and  instead,  works  through  the  interchange  of  wet-surface 
variables  only  that  are  passed  back  and  forth  between  the  DAA  equation  solver  and  that  of 
the  structural  analyzer. 

The  DAA  approach  to  underwater  shock  analysis  involves  the  exact  fluid-structure 
interaction  relations  in  the  limit  of  zero  and  infinite  frequency  while  making  a  smooth  transition 
between  these  limits  in  the  mid-frequency  range.  The  most  typically  used  members  of  the  DAA 
family  available  to  the  user  are  the  one  first  order  DAAl,  and  two  variants  of  a  second  order 
more  accurate  DAA2.  These  latter  are  a  so-called  "modal”  form  DAA2m,  and  a  "curvature 
corrected”  form  DAA2c.  Both  the  DAAl  and  the  DAA2m  use  symmetric  matrices  in  their 
formulations  while  the  DAA2c  requires  a  nonsymmetric  solver  package.  In  terms  of  increasing 
accuracy,  computation  time,  and  memory  requirements  these  choices  rank  in  the  order  DAAl, 
DAA2m,  and  DAA2c,  respectively.  For  special  purposes,  the  limiting  forms  of  the  DAAl, 
the  high  frequency  Plane  Wave  Approximation  (PWA),  and  the  low  frequency  Virtual  Mass 
Approximation  (VMA)  are  also  available  to  the  user. 

A  preprocessor  of  the  USA  code,  FLUMAS,  constructs  a  fluid  mass  matrix  that  fully 
couples  all  of  the  fluid  node  points  on  the  wet  surface  of  the  structure  to  describe  the  low 
frequency  asymptote  of  the  DAA  through  Newton’s  law  of  motion.  This  matrix  is  symmetric 
and  positive  definite  and  is  produced  by  a  boundary  element  formulation  of  Laplace’s  equa¬ 
tion  for  the  irrotational  flow  generated  in  an  inviscid,  incompressible  fluid  by  motions  of  the 
structure.  The  fluid  mass  matrix  also  includes  the  primary  eflfects  of  element  curvature  that 
enhance  the  accuracy  of  the  low  frequency  response. 

The  USA  code  incorporates  a  variety  of  symmetry  and  boundary  conditions  in  its 
wet-surface  fluid  model.  Half  and  quarter  models  are  easily  treated  to  save  computational 
resources  for  symmetric  problems.  In  addition,  free  surface  effects  are  included  through  an 
imaging  process.  The  code  also  has  a  bottom  effects  model  but  this  can  only  be  regarded  as 
a  crude  approximation  to  reality  in  the  absence  of  any  real  information. 
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9.3  Structual  Response  Equation 

Xli6  discrstizBci  diffisrsiiticil  equation  system  for  the  dynamic  response  of  a  linear  structure  can 
be  expressed  in  the  form  [10] 


MgX  +  CsX  +  KgX  —  f 

where  a;  is  a  structural  displacement  vector,  Mg,  C  and  Kg  are  the  symmetric  linear 
structural  mass,  damping  and  stiffness  matrices,  respectively.  /  is  the  external  force  vector, 
and  a  (■)  denotes  temporal  derivative. 

For  excitation  of  a  submerged  structure  by  an  acoustic  wave,  /  is  given  by 


/  =  -GAf{pj+Pg)  +  fD 

where  pj  and  Pg  are  nodal  pressure  vectors  for  the  wet-surface  fluid  mesh  pertaining 
to  the  (known)  incident  wave  and  the  (unknown)  scattered  wave,  respectively,  /o  is  the  dry 
stucture  applied  force  vector,  Af  is  the  diagonal  area  matrix  associated  with  elements  in  the 
fluid  mesh,  and  G  is  the  transformation  matrix  that  relates  the  structural  and  fluid  surface 
forces. 

The  Doubly  Asymptotic  Approximation  may  be  written  [50,51] 


M  fPg  +  pcAfPs  =  pcMfiig 

where  Ug  is  the  vector  of  scattered-wave  fluid  particle  velocities  normal  to  the  structure’s 
wet  surface,  p  and  c  are  the  density  and  sound  velocity  of  the  fluid,  respectively,  and  M/  is 
the  symmetric  fluid  mass  matrix  for  the  wet-surface  fluid  mesh. 

For  excitation  by  an  incident  acoustic  wave,  Ug  is  related  to  structural  response  by  the 
kinematic  compatibility  relation 


G^x  =  w/  -h  «s 


Generally,  G  is  a  rectangular  matrix  whose  height  greatly  exceeds  its  width,  inasmuch 
as  the  number  of  structural  DOF  usually  considerably  exceeds  the  number  of  fluid  DOF.  G 
is  constructed  such  that  only  the  translational  DOF  for  the  structure  couple  with  the  fluid 

DOF. 

Combining  the  above  equations  yields  the  interaction  equations 
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MsX  +  CsX  +  KsX  =  -GAf{pj  +  pj 


M fPg  +  pcAfPg  —  pcM /{G'^x  -  ui) 

These  equations  may  be  solved  simultaneously  at  each  step,  but  such  a  procedure  is 
exceedingly  expensive,  because  of  the  large  connectivity  of  the  coefScient  matrix  involved. 
Efficient  computation  is  possible  trough  the  application  of  a  staggered  solution  procedure, 
that  is  em  unconditionally  stable  with  respect  to  the  choice  of  time  increment,  at  least  for 
linear  structural  response.  Staggered  solution  procedure  is  effected  from  the  above  equations 
(assuming  that  Ms  is  nonsingular)  by  extracting  G^x  from  the  first  equation  into  the  latter 
one.  Premultiplication  of  the  resulting  equation  by  AfM f~^  then  yields 


Afp,  +  {Dfi+  Ds)Ps  =  -pcAfG'^M,  ^{CsX  +  K^x)  -  DsPj  -  pcAfUj 


where 


Dfi  =  pcAfMf  ^Af 


Ds  =  pcAfG'^Ms~^GAf 

Note  that  both  D/i  and  D*  are  symmetric.  This  process  of  injecting  one  of  the  coupled 
equations  into  the  other  to  achieve  stability  is  termed  augmentation,  hence  the  above  is  called 
are  herein  termed  augmented  interaction  equations. 

The  solution  of  the  above  equation  is  carried  in  the  USA  routines,  while  PHLEX  code 
provides  time  integration  of  the  elastic  structure  (with  the  crack)  and  provides  to  the  USA 
module  the  value  of  the  vector  {CsX  +  K^x)  with  possible  dry  structure  force  fj^.  USA 
module  in  turn  returns  the  vector  of  wet  surface  pressures  {pj  +Pg) 

9.4  USA-PHLEXcrack  Interface 

At  this  time  the  interface  between  USA  and  PHLEXcrack  incorporates  only  the  capability  to 
treat  structural  models  that  are  constructed  with  quad  faces  on  the  wet-surface  boundary.  Tri¬ 
angular  faces  are  treated  as  a  subset  of  the  quad  data  structures.  The  Surface- Of-Revolution 
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(SOR)  simplification  for  structural  beam  models  has  not  been  implemented  nor  has  the  ca¬ 
pability  to  treat  cavitation  in  the  surrounding  fluid.  Coupling  of  a  general  DAA  boundary 
to  a  beam  model  of  a  structure  for  whipping  analysis  has  also  been  left  out  of  the  current 
interface.  This  initial  effort  to  couple  USA  with  PHLEXcrack  has  concentrated  on  only  the 
most  basic  type  of  analysis  to  demonstrate  the  utility  of  the  combined  software  package.  The 
above  features  of  the  USA  code  can  be  added  to  the  interface  package  as  the  need  arises. 

The  interface  allows  for  different  time  steps  to  be  used  in  USA  and  PHLEXcrack 
with  all  necessary  interpolations  carried  out  on  the  USA  side  of  the  interface.  Typically  the 
PHLEXcrack  time  step  may  be  smaller  than  the  USA  step  size. 

Special  class  of  boundary  conditions,  called  “WetFace,”  has  been  created  in  PHLEX¬ 
crack  for  this  purpose.  Those  EC’s  should  be  defined  on  element  faces  defining  the  surface  of 
the  body  that  will  be  subjected  to  the  shock  wave.  The  amount  of  pressure  applied  on  these 
surfaces  are  determined  interactively  from  USA  code. 

PHLEXcrack  is  dimensionally  neutral,  i.e.  it  accepts  quantities  and  mesh  dimensions 
in  any  unit  system,  as  long  as  all  of  them  are  in  the  defined  in  the  consistent  way.  USA  code, 
on  the  other  hand,  being  CFD  code,  requires  that  the  quantities  are  defined  in  dimensionless, 
normalized  units.  In  the  examples  presented  in  the  manual  we  assumed  the  characteristic  size 
of  the  model  (outer  radius  of  the  shell),  fluid  density,  and  shock  wave  speed  in  the  fluid  as  all 
equal  1.0.  Thus  the  unit  of  time  is  half  of  the  amount  of  time  it  takes  for  the  shock  wave  to 
travel  in  fluid  from  top  to  the  bottom  of  the  cylinder,  the  density  of  the  shell  material  is  7.85, 
etc. 


9.5  Computational  Procedure 

Usage  of  the  combined  USA-PHLEXcrack  software  involves  the  following  multi-step  procedure: 

1. )  Preparation  of  a  structural  model 

2. )  Run  PHLEXcrack  in  a  pre-processing  mode  to  create  a  database  file  that  is  accessed 
by  two  USA  preprocessors,  FLUMAS  and  AUGMAT 

3. )  Run  FLUMAS  to  create  a  fluid  mass  matrix  appropriate  to  the  wet-surface  geometry 
of  the  structural  model  as  well  as  a  variety  of  subsidiary  fluid  data 

4. )  Run  AUGMAT  to  create  the  so-called  ’’augmented  matrices”  that  are  used  in  the 
computational  form  of  the  DAA 

5. )  Restart  the  USA-PHLEXcrack  processor  to  conduct  the  underwater  shock  analysis 
with  incident  wave  data  and  explosive  charge  location  supplied  by  USA  input  records 

6. )  Post-process  transient  response  results  using  PHLEXcrack  utilities  and  the  USA 
post-processor  POSTER 
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Input  records  for  all  USA  processors  are  listed  in  a  series  of  manuals  that  describe 
each  variable  and  the  order  in  which  they  must  appear  in  the  input  files.  The  POSTER 
utility  incorporates  the  capability  to  compute  pseudo- velocity  shock  spectra  that  is  of  interest 
to  those  who  design  equipment  mounts  and  is  not  generally  found  in  a  typical  off-the-shelf 
post-processors. 


9.6  Software  Validation 

Two  standard  benchmark  problems  have  been  studied  to  test  the  new  interface.  These  are  the 
infinite  cylindrical  shell  and  the  spherical  shell,  both  excited  by  a  plane  step  wave.  Closed  form 
modal  solutions  for  both  problems  have  been  discussed  in  the  literature  and  the  USA  code 
has  demonstrated  its  capabilities  by  replicating  those  solutions  with  a  variety  of  linear  and 
nonlinear  structural  codes  in  the  past.  Computations  with  the  USA-PHLEXcrack  software 
exhibit  the  same  excellent  comparisons  with  these  closed  form  solutions  and  other  computa¬ 
tional  simulations  with  only  the  typical  differences  caused  by  a  variety  of  element  formulations 
used  by  each  structural  code. 


10  Mixed  Elements  in  PHLEXcrack:  EPM  and  Classi¬ 
cal  Shell  Elements 


PHLEXcrack  can  handle  both  EPM  three  dimensional  elements  and  the  shell  elements  com¬ 
bined  together  in  the  same  model.  Being  meshless,  EPM  elements  are  naturally  suited  for 
modeling  crack  propagation.  The  advantages  are  that  cracks  can  advance  in  any  direction  and 
are  not  confined  to  element  faces.  While  shell  elements  are  reserved  to  model  areas  of  thin 
plates.  Both  static  and  dynamic  analyses  can  be  conducted  on  the  combined  model. 


11  User  Routines  for  Crack  Propagation 

The  computer  program  delivered  for  this  project  is  setup  in  such  a  way  to  allow  the  user  to 
add  his  or  her  own  routines  for  crack  propagation.  The  user  can  add  up  to  two  routines  for 
the  calculation  of  fracture  quantities  (SIF’s)  and  up  to  two  routines  for  propagation  physics. 

The  outline  of  these  routines  has  been  coded  and  is  delivered  along  with  the  code. 
In  addition,  the  source  code  for  routines  of  the  fracture  methods  and  physics  described  in 
Sections  5  and  6  are  delivered.  They  are  intended  to  be  used  for  guidance  in  writing  routines 
for  similar  purposes. 

Once  a  new  routine  is  written,  it  has  to  be  compiled  and  linked  with  the  code.  The 
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user  can  then  use  the  new  routine  by  setting  some  parameters  in  the  tcl  file.  Please  refer  to 
the  User  Manual  for  information  on  these  parameters. 


12  Controlling  the  accuracy  of  numerical  algorithms 


Based  on  limited  set  of  numerical  experiments  (some  of  them  are  presented  in  section  13)  we 
present  here  some  practical  notes  on  the  limitations  of  the  code,  and  on  proper  selection  of 
various  control  parameters.  In  particular  element  sizes,  although  significantly  larger  than  in 
traditional  FEM  models,  should  be  selected  appropriately  to  local  crack  accuracy. 


12.1  Crack  details  and  element  sizes 

Although  EPM  allows  to  place  crack  arbitrarily  inside  the  element,  and  it  calculates  proper 
element  stiffness  matrices  (or  their  equivalents  in  EPM),  it  cannot  properly  account  for  all 
minute  details  of  crack  shape  within  element.  Our  experiments  suggest  that  any  details  of 
crack  surface  (sudden  turns,  sharp  comers,  etc.)  which  are  smaller  than  approximately  1/10 
of  the  element  size  are  neglected  in  the  solution.  Crack  curvature  should  not  be  larger  than 
comparable  to  element  size,  as  well  as  crack  having  large  variations  of  curvatures  within  single 
element,  will  be  not  solved  with  very  high  accuracy. 


12.2  Crack  shape  and  accuracy  of  SIF  estimation 

As  explained  in  section  5,  stress  intensity  factors  are  computed  on  the  basis  of  averaged  stress 
in  the  vicinity  of  crack  front,  and  are  based  on  the  assumption  that  the  crack  within  this 
vicinity  will  be  closely  approximated  by  the  planar,  two-dimensional  problem.  The  size  of  the 
domain  for  the  SIF  calculation  (see  parametersdescription  in  “PHLEXcrack  manual”)  should 
be  adjusted  appropriately  to  local  curvatures  of  the  crack  surface.  The  length  of  the  cylindrical 
domain  is  adjusted  automatically  (to  the  length  of  the  straight  segment  of  crack  front),  but 
the  diameter  of  the  cylinder  and  other  parameters  have  to  be  selected  by  the  user. 


12.3  Crack  smoothness 

In  order  to  improve  stability  of  crack  growth,  and  limit  the  restrictions  presented  in  the  last 
section,  PHLEXcrack  contains  several  algorithms  ensuring  proper  smoothness  of  the  crack: 

•  the  angle  between  to  consecutive  crack  velocity  vectors  is  limited. 
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•  velocities  between  neighboring  points  along  crack  edge  are  smoothed  out  using  weighted 
averaging  between  neighbors, 

•  lengths  of  the  segments  along  the  crack  front  are  adjusted  to  fit  between  user  prescribed 
bounds,  and  also  based  on  the  local  curvature  of  the  crack  front. 

To  see  results  of  automatic  “refinement”  of  the  crack  front,  see  bottom  and  right  parts  of  the 
crack  presented  in  the  Fig.  54 


12.4  Crack  folds 

Crack  propagation  algorithm  is  similar  to  the  explicit  time  integration  method  (next  position 
of  the  crack  front  is  calculated  based  on  the  stresses  computed  at  current  crack  position). 
Fortunately,  time  stepping  required  for  the  dynamic  fracture  is  usually  sufficiently  small,  and 
the  problem  may  arise  in  pseudo-static  analysis  (i.e.  without  accurate  modeling  of  shock 
waves).  There  is  however  another  similar  restriction,  resulting  from  the  discrete  (segmented) 
approximation  of  the  crack  front.  In  general,  user  should  adjust  parameters  to  ensure  that 
crack  movement  within  one  time  step  (velocity  times  delta  time)  is  smaller  than  the  minimal 
crack  segment  length.  Fig.  54,  top  portion,  shows  exactly  this  case,  when  crack  propagation 
vectors  were  too  long  compared  with  segment  length  and  crack  “folded”  on  the  small  scale 
(comparable  with  the  segment  length).  PHLEXcrack  geometry  engine  is  not  designed  to 
handle  such  geometries,  and  this  type  of  results  will  most  probably  cause  internal  code  error 
and  crash. 

Crack  engine  can  handle  properly  intersection  of  the  crack  with  the  boundary,  or  with 
other  parts  of  the  crack.  When  such  intersection  is  detected,  the  crack  front  is  stopped  at  the 
node,  split  into  smaller  multiple  fronts  if  necessary,  and  solution  can  proceed.  User  should 
however  verify,  that  the  above  described  accuracy  and  smoothness  requirements  are  satisfied, 
otherwise  further  computations  may  have  no  physical  meaning. 


13  Numerical  Examples 

The  cracked  element  partition  method  (CEPM)  described  previously  is  used  in  this  section  to 
solve  several  illustrative  examples.  In  all  examples,  the  stress  intensity  factors  are  computed 
using  the  technique  outlined  in  Section  5. 
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13.1  Single  Crack  with  Mode  I  Solution  Under  Static  Loading 

In  this  section,  the  edge-cracked  panel  illustrated  in  Fig.  11  is  analyzed  using  the  CEPM.  In 
the  computations,  the  following  parameters  are  assumed:  h  —  b  —  1.0,  a  =  0.5,  distributed 
tractions  a  =  1.0  and  uniform  thickness  t  =  0.1.  The  material  is  assumed  to  be  linearly  elastic 
with  E  =  1000.0  and  u  =  0.3. 


Figure  11:  Single  edge-noth  test  specimen 

The  domain  is  discretized  using  the  hexahedral  mesh  shown  in  Fig.  12.  There  are  961 
elements  in  the  mesh.  A  state  of  plane  strain  is  modeled  by  constraining  the  displacement  in 
the  z-direction  at  2  =  0  and  z  =  t.  The  representation  of  the  crack  surface  is  shown  in  Fig.  13. 
It  is  composed  of  four  triangles  and  two  edge  elements  (used  to  identify  the  crack  front).  Note 
that  these  triangles  are  used  only  for  the  geometric  representation  of  the  crack  surface  and  that 
there  is  no  degrees  of  freedom  associated  with  then.  The  stress  intensity  factors  reported  for 
this  problem  are  computed  at  aj  =  (0.5, 1.0, 0.05)  which  are  the  coordinates  of  a  node  of  an  edge 
element  and  is  located  at  the  middle  plane  of  the  body.  The  geometric  definition  of  the  outer 
skin  of  the  domain  (shown  on  Fig.  13)  is  also  given  as  input  data  for  the  geometric  engine. 
The  base  vectors  of  the  coordinate  system  associated  with  singular  functions  used  at  the  crack 
front  is  also  displayed  in  Fig.  13.  The  base  vectors  and  corresponding  singular  functions  are 
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computed  completely  automatically  using  the  geometric  engine.  This  functionally  is  specially 
important  during  dynamic  crack  propagation  simulations  or  when  the  geometry  of  the  domain 
or  crack  surface  are  not  so  trivial  as  in  this  example. 

Figure  14  shows  a  closer  look  at  the  discretization  near  the  crack  front.  It  can  be 
observed  that  the  crack  surface  does  not  respect  the  element  boundaries-/^  can  arbitrarily  cut 
the  elements  in  the  mesh.  The  nodes  carrying  singular  degrees  of  freedom  are  represented 
by  diamond-shaped  dots.  The  singular  functions  used  at  these  nodes  are  those  presented  at 
Sections.  Whenever  they  are  used,  only  the  nodes  of  the  elements  that  contain  the  crack  front 
are  enriched  with  these  singular  shape  functions  (in  this  example,  there  is  only  one  element  at 
the  crack  front).  The  enrichment  of  the  elements  at  the  crack  front  with  appropriate  singular 
functions  is  done  automatically  using  the  geometric  engine  to  construct  appropriate  coordinate 
systems  at  the  crack  front. 

The  computed  values  oi  Kj  and  of  the  strain  energy,  U,  for  several  discretization  are 
shown  on  Tables  1  and  2.  In  the  tables,  N  denotes  the  number  of  degree  of  freedom  of 
a  particular  discretization.  The  discontinuity  in  the  displacement  field  for  those  elements 
completely  severed  by  the  crack  surface  is  modeled  using  the  technique  described  in  Section 
2.4.  For  the  element  at  the  crack  front,  two  approaches  are  used.  In  the  first  case,  the  visibility 
approach  in  combination  or  not  with  singular  functions  at  the  crack  front  is  used.  The  results 
using  this  approach  is  shown  in  Table  1.  In  the  second  case,  the  wrap-around  approach  in 
combination  with  singular  functions  at  the  crack  front  is  used.  The  results  using  this  approach 
is  shown  in  Table  2.  The  notation  p  =  H-p^,  1  +Py,  1  +^2  =  1  +  iPx-,Py,Vz)  is  used  to  denote 
the  polynomial  order  of  the  approximation  over  non-cracked  elements.  The  “I’s”  indicating 
the  linear  order  of  the  functions  defining  the  partition  of  unity  (linear  hexahedral  finite  element 
shape  functions  in  this  case)  and  Pa,,  Py,Pz  denote  the  degrees  of  the  polynomial  basis  functions 
Lia  in  the  x,y,z  directions,  respectively  (see  Section  2.3  for  the  definition  of  Lia).  Therefore, 
if  only  the  partition  of  unity  is  used,  we  have  p  =  l-|-0, 1-f-0, 1-|-O  =  l  +  0.  A  quadratic 
approximation  in  the  plane  XY  and  linear  in  the  ^-direction  is  denoted  byp  =  H-l,l-l-l,l-l- 
0  =  14-  (1,1,0).  For  the  cracked  elements,  the  polynomial  degree  of  the  approximation  is 
denoted  by  Pc  =  O+p^,  0+Py,  0+Pz  =  0+{pa;,Py,Pz).  The  “O’s”  indicating  that,  in  general,  the 
functions  defining  the  partition  of  unity  over  cracked  elements  can  only  represent  a  constant 
function.  A  quadratic  approximation  in  the  plane  XY  and  linear  in  the  ^-direction  over  cracked 
elements  (fully  or  partially  cracked)  is  denoted  by  Pc  =  0  4-  2, 0  4-  2, 0  -f  1  =  0  -f  (2, 2, 1).  The 
stress  intensity  factors  are  computed  using  the  technique  outlined  in  Section  5.  The  following 
parameters  are  used  in  all  computations  presented  in  this  section: 


•  Dimensions  of  the  extraction  domain:  d  =  (0.4,1.0,0.05).  Which  represents  a  cylinder 
of  radius  0.4  and  length  0.05. 

•  Number  of  integration  points  in  the  r,  6  and  2:  directions:  n  =  (10, 20, 1),  respectively. 

•  Type  of  D  matrix  equal  identity  matrix 
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Figure  12;  CEPM  discretization  for  the  model  problem  shown  on  Fig.  11. 
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Figure  14:  Zoom  at  the  crack  front. 
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•  Power  of  the  weighting  function:  6. 


As  81  reference,  the  value  of  Kj  computed  by  Xada  at  al.  [48]  using  a  boundary  technique 
applied  to  a  two  dimensional  idealization  of  the  problem  shown  on  Fig.  11.  The  value  = 

3.54259  reported  in  [48]  has  an  error  smaller  than  0.5%. 

The  discretizations  Vis-1,  Vis-2  and  Vis-3  do  not  use  singular  functions  at  the  crack 
front  while  the  discretizations  Vis-2,  Vis-4  and  Vis-6  do.  All  ’Vis’  discretizations  use  the 
visibility  approach  to  model  the  crack  at  partially  cut  elements.  It  can  be  observed  from 
Table  1  that  the  use  of  singular  functions  gives  a  noticeable  improvement  on  the  computed 
stress  intensity  factors  while  the  increase  in  the  number  of  degrees  of  freedom  is  only  marginal 
(less  than  one  percent  in  the  case  of  the  discretization  Vis-5).  In  contrast,  the  p-enrichment 
of  the  cracked  elements  only  (discretizations  Vis-1  and  Vis-3)  gives  little  improvement  on  the 
computed  Ki.  Nonetheless,  the  enrichment  does  improve  the  computed  strain  energy  by  about 
3%.  This  behavior  indicates  that  the  technique  used  to  compute  the  stress  intensity  factor  is 
not  optimal  since  optimally  the  computed  stress  intensity  factors  must  converge  at  the  same 
rate  as  the  computed  strain  energy  [46,47]. 

The  discretizations  WA-1,  WA-2  and  WA-3  use  singular  functions  in  combination  Avith 
the  wrap-around  technique  to  model  the  the  crack  at  partially  cut  elements.  Comparing  the 
results  for  the  discretizations  Vis-2  with  WA-1  or  Vis-4  with  WA-2,  it  can  be  observed  that, 
for  the  same  number  of  degrees  of  freedom,  the  wrap-around  approach  gives  better  results 
for  the  stress  intensity  factors  than  the  visibility  approach.  This  is  in  spite  of  the  fact  that 
the  discretizations  using  wrap-around  has  a  smaller  strain  energy  than  the  corresponding 
discretizations  with  using  visibility.  This  can  be  explained  by  the  fact  that  the  visibility 
approach  creates  spurious  discontinuities  in  the  displacement  field  near  the  crack  front  which 
results  in  a  less  stiff  discretization  compared  with  the  wrap-around  approach  (which  does  not 
creates  such  spurious  discontinuities).  It  can  be  observed  that  the  discretization  WA-2  gives 
a  better  value  for  Kj  than  the  discretization  WA-3  in  spite  of  the  fact  that  the  later  giver 
a  larger  value  for  the  strain  energy  than  the  latter.  This,  again,  points  to  limitations  of  the 
technique  used  to  compute  the  stress  intensity  factor. 

Figure  15  shows  a  contour  plot  of  the  displacement  in  the  vertical  direction  near  the 
crack  computed  using  the  discretization  WA-3.  The  discontinuity  in  the  displacement  field 
constructed  using  the  technique  presented  in  Section  2.4  is  clearly  observed.  Figure  16  shows 
a  contour  plot  for  the  von  Vises  stress  computed  with  the  same  discretization  and  Fig.  17 
shows  a  closer  look  near  the  crack  front.  The  computed  stresses  are  all  raw  stresses  computed 
at  arbitrary  points  inside  each  element.  Figure  18  shows  the  same  quantity  computed  using 
the  discretization  vis-5.  It  can  be  observed  that  the  stress  field  is  quite  disturbed  near  the 
crack  front.  This  is  caused  by  the  spurious  discontinuities  created  by  the  visibility  approach 
near  the  crack  front. 
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Table  1:  CEPM  using  visibility  for  the  elements  at  the  crack  front.  The  stress  intensity  factor 
is  computed  at  (0.5,1.0,0.05). 


Discret. 

p  =  1-f- 

Pc  —  0+ 

Sing  Fn 

N 

U  X  10^ 

Ki 

Vis-1 

(0,0,0) 

(1.1,1) 

No 

6,756 

2.27136 

3.2387 

0.91422 

Vis-2 

(0,0,0) 

(1.1.1) 

Yes 

6,900 

2.32941 

3.3741 

0.95244 

Vis-3 

(0,0,0) 

(2,2,1) 

No 

7,776 

2.33848 

3.2495 

0.91727 

Vis-4 

(0,0,0) 

(2,2,1) 

Yes 

7,920 

2.36701 

3.3938 

0.95800 

Vis-5 

(1,1,0) 

(2,2,1) 

No 

19,656 

2.38643 

3.4339 

0.96932 

Vis-6 

(1,1,0) 

(2,2,1) 

Yes 

19,800 

2.42281 

3.4584 

0.97623 

Table  2:  CEPM  using  wrap-around  for  the  elements  at  the  crack  front.  The  stress  intensity 
factor  is  computed  at  (0.5,1.0,0.05). 


Discret. 

+ 

rH 

It 

II 

o 

H- 

i _ 

Sing  Fn 

N 

U  X  10^ 

Ki 

WA-1 

(0,0,0) 

(1.1.1) 

Yes 

6,900 

2.22465 

3.3886 

0.95653 

WA-2 

(0,0,0) 

(2,2,1) 

Yes 

7,920 

2.25596 

3.4228 

0.96618 

WA-3 

(1,1,0) 

(2,2,1) 

Yes 

19,800 

2.29396 

3.3431 

0.94369 

Figure  15:  Displacement  in  the  ?/-direction  computed  with  discretization  WA-3. 
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Figure  16:  Von  Mises  stress  computed  with  discretization  WA-3. 


Figure  17:  Zoom  at  the  crack  front  showing  von  Mises  stress  computed  with  discretization 
WA-3. 
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Figure  18:  Zoom  at  the  crack  front  showing  von  Mises  stress  computed  with  discretization 
vis-5. 
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13.2  An  Inclined  Crack  Problem 


As  another  test  problem,  we  consider  the  cracked  panel  shown  in  Fig.  19.  This  problem  has 
been  analyzed  by  Szabd  and  Babuska  [46]  using  the  p-version  of  the  finite  element  method 
and  by  Oden  and  Duarte  [37]  using  the  hp  Cloud  method.  In  both  references,  plane  stress 
condition  and  unity  thickness  are  used.  Here,  the  plane  stress  condition  is  approximated  by 
using  a  small  thickness,  t  =  0.1,  for  the  domain  compared  to  other  dimensions.  In  addition, 
we  assume  Young’s  modulus  E  =  l,  Poisson’s  ratio  v  =  0.3,  distributed  traction  a  =  1.0  and 
w  =  1  (see  Fig.  19).  These  same  values  are  used  in  references  [46]  and  [37].  We  adopt  as 
a  reference,  the  values  of  Ki,  Kn  and  strain  energy,  U,  computed  by  Oden  and  Duarte  [37]. 
They  are,  respectively, 

=  1.508284 
=  -0.729706 

[/He/  ^  0.170402 

These  values  agree  very  well  with  those  computed  by  Szabo  and  Babuska  [46]  (less  than  0.1% 
difference).  The  value  of  the  was  scaled  to  take  into  account  the  difference  in  thickness 
used  here  {t  =  0.1)  and  adopted  by  Oden  and  Duarte  [37]  {t  =  1.0). 

The  discretization  of  the  domain  using  765  hexahedral  elements  is  shown  on  Fig.  20. 
The  inclined  crack  is  also  shown  in  the  figure.  Figure  21  shows  a  closer  look  near  the  crack. 
It  can  be  observed  that  the  crack  surface  cuts  the  elements  in  the  mesh  in  a  quite  arbitrary 
manner.  In  fact,  the  meshing  of  the  domain  is  done  as  if  there  is  no  crack  at  all.  The  only 
consideration  used  during  the  meshing  of  the  domain  was  to  use  a  more  refined  mesh  near  the 
location  of  the  crack  front.  The  crack  representation  is  created  and  passed  to  the  geometric 
engine  as  input  data.  The  geometric  engine  uses  no  information  whatsoever  about  the  mesh. 
Nodes  in  the  mesh  that  are  too  close  to  the  crack  surface  are  then  automatically  moved  a  small 
distance  away  from  the  crack  surface  (this  can  be  observed  in  Fig.  21  near  the  crack  front). 
This  is  required  because  an  approximation  node  must  be  located  at  one  or  another  side  of  the 
crack  surface  (but  not  (numerically)  on  it).  The  representation  of  the  crack  surface  used  here  is 
topologically  identical  to  the  one  used  in  the  previous  example  (see  Fig.  13).  Figure  21  shows 
a  closer  look  at  the  mesh  and  crack  representation  near  the  crack  front.  The  nodes  carrying 
singular  degrees  of  freedom  are  represented  by  diamond-shaped  dots.  The  singular  functions 
used  at  these  nodes  are  those  presented  at  Section3.  As  in  the  previous  example,  whenever 
they  are  used,  only  the  nodes  of  the  elements  that  contain  the  crack  front  are  enriched  with 
these  singular  shape  functions. 

Figure  23  shows  a  contour  plot  of  the  displacement  in  the  vertical  direction  near  the 
crack  computed  using  the  discretization  WA-3.  The  discontinuity  in  the  displacement  field 
constructed  using  the  technique  presented  in  Section  2.4  is  clearly  observed.  Figure  24  shows 
a  contour  plot  for  the  von  Vises  stress  computed  with  the  same  discretization  and  Fig.  25 
shows  a  closer  look  near  the  crack  front.  The  computed  stresses  are  all  raw  stresses  computed 
at  arbitrary  points  inside  each  element.  Figure  26  shows  the  same  quantity  computed  using 
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the  discretization  vis-5.  It  can  be  observed,  as  in  the  previous  example,  that  the  stress  field 
is  quite  disturbed  near  the  crack  front.  This  is  caused  by  the  spurious  discontinuities  created 
by  the  visibility  approach  near  the  crack  front. 

A  summary  of  the  results  is  presented  on  Tables  3,  4,  5  and  6.  The  notation  used  to 
describe  the  various  discretizations  used  is  the  same  as  in  the  previous  section.  The  stress 
intensity  factors  are  computed  at  a  point  in  the  crack  front  located  at  the  middle  surface  of 
the  body.  The  following  parameters  are  used  for  extracting  the  stress  intensity  factors: 

•  Dimensions  of  the  extraction  domain:  d  =  (0.2,1.0,0.05).  Which  represents  a  cylinder 
of  radius  0.4  and  length  0.05. 

•  Number  of  integration  points  in  the  r,d  and  z  directions:  n  =  (10,40, 1),  respectively. 

•  Type  of  D  matrix  equal  identity  matrix 

•  Power  of  the  weighting  function:  6. 
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Figure  20:  CEPM  discretization  for  the  model  problem  shown  on  Fig.  19. 
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Figure  21:  Zoom  showing  the  elements  cut  by  the  crack  surface. 


Figure  22:  Three-dimensional  view  of  the  crack  front. 
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Figure  23:  Contour  plot  of  the  displacement  in  the  vertical  direction  computed  using  the 
discretization  WA-3. 


Mixed  mode  crack 


Color:  Von  Mises  str. 

Figure  24:  Contour  plot  for  the  von  Vises  stress  field  computed  with  discretization  WA-3. 
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Figure  25:  Zoom  showing  the  von  Mipes  stress  near  the  crack  front  computed  with  discretiza¬ 
tion  WA-3. 


Figure  26:  Zoom  showing  the  von  Mises  stress  near  the  crack  front  computed  with  discretiza¬ 
tion  vis- 5. 


Table  3  shows  the  computed  strain  energy  for  the  various  discretizations  using  the 
visibility  approach  at  the  crack  front  with  or  without  singular  functions.  It  can  be  observed 
that  the  p-enrichment  of  the  approximation  has  a  more  significant  effect  on  the  strain  energy 
values  than  the  addition  of  singular  functions  at  the  crack  front.  Nonetheless,  as  in  the  previous 
example,  the  addition  of  singular  functions  improves  considerably  the  computed  stress  intensity 
factors  (Cf.  Table  4).  In  the  case  of  the  discretization  Vis-3,  for  example,  the  enrichment  with 
singular  functions  adds  only  1.8%  more  degrees  of  freedom  while  the  error  on  the  compute 
value  of  Ki  decreases  from  14.0%  to  only  3.5%  and  the  error  on  the  computed  Kn  decreases 
from  10.8%  to  only  1.4%.  That  is,  the  error  on  the  computed  Kj  and  Kn  decrease  by  75.0% 
and  87,0%,  respectively. 

The  results  for  the  discretizations  that  use  the  wrap-around  approach  at  the  crack 
front  are  presented  in  Tables  5  and  6.  The  discretization  WA-3  has  a  relative  error  in  energy, 
(f/fle/  _  U)/U^^,  of  only  0.02%  which  corresponds  to  a  relative  error  in  the  energy  norm  of 
only  1.41%.  This  same  problem  was  also  solved  using  the  classical  hp  finite  element  method 
with  the  hp  adaptation  driven  by  error  indicators  based  on  the  element  residual  method.  The 
results  obtained  after  seven  adaptive  cycles  are  shown  on  Table  7.  The  relative  error  in  energy 
and  in  the  energy  norm  for  this  discretization  are  0.3%  and  5.5%,  respectively.  Note  that 
this  discretization  has  more  degrees  of  freedom  (18585)  than  the  discretization  WA-3  (16632) 
but  an  error  in  the  energy  norm  almost  four  times  bigger.  The  reason  for  this  is  that  the 
CEPM  discretization  can  capture  the  singular  field  near  the  crack  front  more  efficiently  by 
using  customized  singular  functions. 

Table  3:  CEPM  using  visibility  for  the  elements  at  the  crack  front.  In  the  table,  “Discr” 
stands  for  discretization,  “Sig  Fn”  stands  for  singular  functions  at  the  crack  front,  “N”  stands 
for  number  of  equations  and  “U”  stands  for  strain  energy. 


Discr 

p  =  1+ 

II 

o 

+ 

Sig  Fn 

N 

U 

[//t/iie/ 

Vis-1 

(0,0,0) 

(1,1.1) 

No 

5874 

0.157986 

0.9271 

Vis-2 

(0,0,0) 

(1,1.1) 

Yes 

6018 

0.158644 

0.9310 

Vis-3 

(0,0,0) 

(2,2,1) 

No 

7764 

0.164454 

0.9651 

Vis-4 

(0,0,0) 

(2,2,1) 

Yes 

7908 

0.165053 

0.9686 

Vis-5 

(1,1,0) 

(2,2,1) 

No 

16488 

0.171148 

1.0044 

Vis-6 

(1,1,0) 

(2.2,1) 

Yes 

16632 

0.171558 

1.0068 
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Table  4:  CEPM  using  visibility  for  the  elements  at  the  crack  front.  In  the  table,  “Discr” 
stands  for  discretization,  “Sig  Fn”  stands  for  singular  functions  at  the  crack  front  and  “N” 
stands  for  number  of  equations. 


Discr 

p  =  1+ 

o 

II 

O 

+ 

SigFn 

N 

Ki 

Kii 

Ki/K^^ 

Vis-1 

(0,0,0) 

(1,1.1) 

No 

5874 

1.1099 

-0.58755 

0.7286 

0.8052 

Vis-2 

(0,0,0) 

(1.1.1) 

Yes 

6018 

1.2656 

-0.68936 

0.8391 

0.9447 

Vis-3 

(0.0,0) 

(2.2,1) 

No 

7764 

1.2964 

-0.65095 

0.8595 

0.8921 

Vis-4 

(0,0,0) 

(2,2,1) 

Yes 

7908 

1.4550 

-0.71979 

0.9647 

0.9864 

Vis-5 

(1,1,0) 

(2,2,1) 

No 

16488 

1.4904 

-0.73524 

0.9881 

1.0076 

Vis-6 

(1,1.0) 

(2,2,1) 

Yes 

16632 

1.6595 

-0.77266 

1.1002 

1.0589 

Table  5:  CEPM  using  wrap-around  for  the  elements  at  the  crack  front.In  the  table,  “Discr” 
stands  for  discretization,  “Sig  Fn”  stands  for  singular  functions  at  the  crack  front,  “N”  stands 
for  number  of  equations  and  “U”  stands  for  strain  energy. 


Discr 

p  =  l+ 

Pc  —  0+ 

Sig  Fn 

N 

U 

u/u^r 

WA-1 

(0,0,0) 

(1,1,1) 

Yes 

6018 

0.157978 

0.9271 

WA-2 

(0,0,0) 

(2,2,1) 

Yes 

7908 

0.164183 

0.9635 

WA-3 

(1,1,0) 

(2,2,1) 

Yes 

16632 

0.170372 

0.9998 

Table  6:  CEPM  using  wrap-around  for  the  elements  at  the  crack  front.In  the  table,  “Discr” 
stands  for  discretization,  “Sig  Fn”  stands  for  singular  functions  at  the  crack  front  and  “N” 
stands  for  number  of  equations. 


Discr 

p  =  1+ 

Pc  =  0+ 

Sig  Fn 

N 

Ki 

Kii 

Ki/Kf^ 

KiilKf^^ 

WA-1 

(0,0,0) 

(1,1.1) 

Yes 

6018 

1.3175 

-0.64084 

0.8735 

0.8782 

WA-2 

(0,0,0) 

(2,2,1) 

Yes 

7908 

1.4663 

-0.69115 

0.9722 

0.9472 

WA-3 

(1,1,0) 

(2,2,1) 

Yes 

16632 

1.6290 

-0.72814 

1.0800 

0.9979 

Table  7:  Results  using  hp  finite  element  method  and  seven  adaptive  cycles. 


Discr 

N 

U 

u/u^^f 

Ki 

Kii 

Ki/Kf^ 

hp  FEM 

18585 

0.16988 

0.9969 

1.4680 

-0.7041 

0.9733 

0.9649 
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13.3  Plate  Under  Impact  Load 


In  this  section  we  investigate  the  performance  of  the  CEPM  in  modeling  propagating  cracks 
in  a  body  subjected  to  impact  loads.  The  test  problem  is  illustrated  on  Fig.  27.  This  problem 
has  been  analyzed  by  Lu  et  al.  [32],  Krysl  and  Belytschko  [24],  Organ  [44]  and  Belytschko  and 
Tabbara  [5]  using  the  element  free  Galerkin  method,  by  Gallego  and  Dominguez  [21]  using  a 
boundary  element  method,  among  others.  A  state  of  plane  strain  and  the  following  parameters 
are  adopted 

•  Dimensions:  b  =  10.0,  h  =  2.0,  a  =  5.0  and  uniform  thickness  t  =  0.1, 

•  Loading:  <T(t)  ^  a  H{t)  -  63750.0  H{t),  t  >  0.  Here,  H{t)  is  the  Heaviside  step 
function. 

•  Material  Properties:  Linear  elastic  material  with  E  —  2.0  x  10^\  u  =  0.3  and  p  —  7833.0. 

•  Time  Step:  At  =  10”®. 

A  state  of  plane  strain  is  modeled  by  constraining  the  displacement  in  the  z-direction 
a,t  z  =  0  and  z  =  t.  Two  uniform  hexahedral  meshes  are  used.  The  first  one  has  125,  49 
and  1  element  in  the  x-,  y-  and  z-directions,  respectively.  This  same  mesh  was  used  in  the 
computations  of  Krysl  and  Belytschko  [24].  We  denote  this  as  the  fine  mesh.  The  second 
mesh  has  65,  25  and  1  element  in  the  x-,  y-  and  z-directions,  respectively  and  is  shown  on 
Fig.  28  along  with  the  crack  representation.  This  mesh  is  denoted  as  the  coarse  mesh.  The 
representation  of  the  crack  surface  and  of  the  outer  skin  of  the  body  are  shown  in  Figs.  29 
and  30.  It  is  composed  of  five  triangles  and  four  edge  elements.  There  are  five  vertex  nodes 
uniformly  spaced  at  the  crack  front.  The  stress  intensity  factors  reported  for  this  problem  are 
computed  at  a;  =  (5.0, 2.0, 0.05)  and  is  located  at  the  middle  plane  of  the  body.  This  vertex 
node  is  denoted  vertex  3  hereafter. 

The  stress  intensity  factors  are  computed  using  the  technique  outlined  in  Section  5. 
The  following  parameters  are  used  in  all  computations  presented  in  this  section: 

•  Dimensions  of  the  extraction  domain:  d  =  (1.0,1.0,0.05).  Which  represents  a  cylinder 
of  radius  1.0  and  length  0.05. 

•  Number  of  integration  points  in  the  r,  6  and  2:  directions:  n  =  (5, 30, 1),  respectively. 

•  Type  of  D  matrix  equal  material  matrix 

•  Power  of  the  weighting  function:  3. 

Four  CEPM  discretizations  are  used  (the  notation  used  for  p  and  Pc  is  defined  in 
Section  13.1)  : 
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•  Discretization  1: 

-  Fine  mesh  (125  x  49  x  1  elements), 

—  Degree  of  approximation  over  non  cracked  elements  p  =  1  +  0. 

—  Degree  of  approximation  over  cracked  elements  Pc  =  0  +  1. 

-  Algorithm  to  build  a  discontinuous  partition  of  unity  over  the  elements  at  the  crack 
front:  Visibility  (and  no  singular  functions). 

•  Discretization  2: 

-  Fine  mesh  (125  x  49  x  1  elements), 

—  Degree  of  approximation  over  non  cracked  elements  p  =  1  +  0. 

—  Degree  of  approximation  over  cracked  elements  pc  =  0  +  1. 

-  Algorithm  to  build  a  discontinuous  partition  of  unity  over  the  elements  at  the  crack 
front:  Wrap-around  and  singular  functions. 

•  Discretization  3: 

—  Coarse  mesh  (65  x  25  x  1  elements), 

—  Degree  of  approximation  over  non  cracked  elements  p  =  1  +  0. 

-  Degree  of  approximation  over  cracked  elements  pc  =  0  +  1. 

-  Algorithm  to  build  a  discontinuous  partition  of  unity  over  the  elements  at  the  crack 
front:  Visibility  (and  no  singular  functions). 

•  Discretization  4: 

-  Coarse  mesh  (65  x  25  x  1  elements), 

—  Degree  of  approximation  over  non  cracked  elements  p  =  1  +  0. 

—  Degree  of  approximation  over  cracked  elements  Pc  =  0 +  (2,2,1). 

-  Algorithm  to  build  a  discontinuous  partition  of  unity  over  the  elements  at  the  crack 
front:  Visibility  (and  no  singular  functions). 


13.3.1  Reference  Solution 

The  CEPM  results  are  compared  to  the  analytic  solution  for  a  semi-infinite  crack  in  the ’plane 
proposed  by  Freund  [20].  The  problem  solved  by  Freund  is  represented  in  Fig.  31.  The 
two-dimensional  domain  has  a  straight  semi-infinite  crack,  is  in  a  state  of  plane  strain  and  is 
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Figure  27:  Model  problem  used  to  analyze  the  performance  of  the  CEPM  in  modeling  propa 
gating  cracks  in  a  body  subjected  to  impact  loads. 


Figure  28:  CEPM  Discretization. 


Figure  30:  CEPM  Discretization. 


loaded  by  uniformly  distributed  tractions  applied  at  time  t  =  0.  The  mode-I  stress  intensity 
factor,  as  a  function  of  time  and  crack  speed  C  is  given  by  (  [20]). 

.. -  ¥(c) . /(I  -  2,/)(i  -  i)c, 

where, 


•  d  is  the  magnitude  of  the  tensile  tractions, 

•  Cd  is  the  pressure  wave  speed  in  the  body  which  is  given  by 

^  _  U(«+i) 

where  /i  =  k  =  Z  -  Au  (for  plane  strain  state).  For  the  material  properties 

given  previously,  we  get  Cd  —  5862.7. 

•  t  is  the  time  the  elastic  wave  hits  the  crack.  For  the  problem  represented  in  Fig.  27  and 
Cd  =  5862.7 

i  =  0.000341 

•  k{C)  is  scaling  factor  that  takes  into  account  that  the  crack  front  is  advancing  with 
speed  C  and  is  given  by 

_  1°  -  C/gH 

'  '“1.0-0.5*C/Cs 

where  Cr  is  the  Rayleigh  wave  speed.  For  this  test  problem,  with  the  material  properties 
given  above,  Cr  =  3030  (see  Section  6). 


Due  to  symmetries,  Kn  =  0.  The  magnitude  of  the  energy  release  rate  is  given  by 

G(C,t)  =  (l-£)G(C  =  0,i) 

where 


and,  for  plane  strain. 


It  should  be  noted  that  due  to  the  finite  dimensions  of  the  domain  modeled  here  there 
will  be  waves  reflected  by  the  boundary.  This  reflected  waves  will  eventually  reach  the  crack 
front  and  a  comparison  of  the  numerical  solution  with  the  above  reference  solution  will  no 
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longer  be  valid.  The  first  refiected  wave  to  reach  the  crack  is  a  pressure  wave  after  traveling 
from  a  loaded  edge  to  the  opposite  edge  and  then  back  to  the  crack  front  [24].  This  happens 

t=  —  =z  — ^  =  0.00102 

Ca  5862.7 

A  more  detailed  discussion  on  the  wave  patterns  that  reach  the  crack  front  can  be  found 
in  [24]. 


13.3.2  Stationary  Crack 

As  a  first  test,  we  consider  the  case  of  a  stationary  crack.  Discretizations  1  and  2  as  described 
above  are  used.  The  computed  mode-I  and  mode-II  stress  intensity  factors  and  the  energy 
release  rate  G  are  plotted  on  Figs.  32  and  33,  respectively.  It  can  be  observed  a  good  agreement 
between  the  computed  and  reference  values  for  t  <  i.  The  finite  dimensions  of  the  extraction 
domain  for  the  stress  intensity  factors  and  of  the  support  of  the  shape  functions  are  responsible 
for  the  non-zero  values  computed  before  the  pressure  wave  hits  the  crack  front  {at  t  =  i).  It 
can  also  be  observed  that  both  discretizations  give  very  close  results. 


13.3.3  Moving  Crack  with  Prescribed  Speed  and  Direction 
Here,  the  crack  speed  is  given  by 

C{t)  =  H{t  -  0.00044)^  =  H{t  -  0.00044)1010.0 

O 
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Time  History,  Plane  Strain  Model 


Figure  32:  Time  history  for  mode-I  and  mode-II  stress  intensity  factors  Kj  and  Kn  using 
Discretizations  1  and  2.  Stationary  crack. 


Time  History,  Plane  Strain  Model 


Figure  33:  Time  history  for  the  energy  release  rate  G  using  Discretizations  1  and  2.  Stationary 
Crack. 
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and  propagates  straightly  along  the  direction  0  =  0  where  0  is  measured  from  the  x-axis 
counter-clockwise.  The  values  of  Kj,  Ku  and  G  computed  with  the  Discretization  1  are 
plotted  in  Figs.  34  and  35  against  the  reference  solution.  While  the  computed  values  present 
some  oscillation,  they  are  in  good  agreement  with  the  reference  solution  before  t  =  t  (when 
reflected  waves  hits  the  crack  surface).  Organ  [44],  Krysl  and  Belytschko  [24]  and  Belytschko 
and  Tabbara  [5]  also  reported  oscillations  on  their  results  obtained  with  the  element  free 
Galerkin  method. 

Time  History,  Plane  Strain  Model 


Figure  34:  Time  history  for  mode-I  and  mode-II  stress  intensity  factors  Ki  and  Ku  using 
Discretization  1.  Moving  Crack  with  prescribed  speed  and  direction. 


13.3.4  Moving  Crack  with  Prescribed  Speed 

Here,  the  direction  of  the  crack  advancement  is  given  by  (10).  Discretization  1  (fine  mesh) 
and  3  (coarse  mesh)  are  used.  Figures  36  and  37  shows  the  time  history  for  Ki,  Ku  and  G, 
respectively.  It  can  be  observed  that  CEPM  solution  was  able  to  capture  very  well  the  slop  of 
the  reference  solution.  Figure  38  shows  the  crack  surface  at  time  t  =  0.00168s. 

The  results  obtained  with  Discretization  3  (coarse  mesh)  are  shown  on  Figs.  39  and 
40.  The  approximation  of  the  slope  of  the  reference  solution  is  still  quite  reasonable,  specially 
while  the  crack  remains  stationary  {t  <  0.00044).  The  results  obtained  with  Discretizations  1 
and  3  are  compared  on  Figs.  41  and  42. 

The  eflfect  of  p  enrichment  of  the  cracked  elements  is  investigated  by  using  Discretization 


Time  History,  Plane  Strain  Model 


Figure  35:  Time  history  for  the  energy  release  rate  G  using  Discretization  1.  Moving  Crack 
with  prescribed  speed  and  direction. 


Time  History,  Plane  Strain  Model 


Figure  36:  Time  history  for  Ki  and  Ku  using  Discretization  1.  Moving  crack  with  advance¬ 
ment  direction  given  by  (10). 
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Time  History,  Plane  Strain  Model 


Figure  37:  Time  history  for  G  using  Discretization  1.  Moving  crack  with  advancement  direc¬ 
tion  given  by  (10). 

4  (which  is  identical  to  Discretization  3  but  uses  pc  =  0-|-(2,2,l)  instead  of  Pc  =  0  +  1.  The 
results  for  this  discretization  are  shown  on  Figs.  43  and  44.  It  can  be  observed  the  the  p 
enrichment  of  the  cracked  elements  increases  the  amplitude  of  the  oscillations  of  the  computed 
quantities.  Figure  45  is  identical  to  Fig.  44  but  here  the  tic  marks  of  the  x-axis  are  placed 
exactly  at  the  times  when  the  crack  front  crossed  the  boundary  between  two  elements.  It  can 
be  observed  that  the  peaks  in  the  oscillations  occur  just  before  those  instants.  Note  that  as 
the  crack  crosses  the  boundary  between  two  elements  it  passes  close  to  the  nodes.  This  same 
phenomena  was  observed  by  Organ  [44]. 
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Figure  38:  Crack  surface  at  t  =  0.00168s.  Moving  crack  with  advancement  direction  given  by 


Time  (sec) 

Figure  39:  Time  history  for  Ki  and  Ku  using  Discretization  3.  Moving  crack  with  advance¬ 
ment  direction  given  by  (10). 


Time  History,  Plane  Strain  Model 


Figure  40:  Time  history  for  G  using  Discretization  3.  Moving  crack  with  advancement  direc¬ 
tion  given  by  (10). 
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Figure  41:  Comparison  of  results  for  Kj  and  Ku  computed  with  Discretizations  1  and  3. 


Time  History,  Plane  Strain  Model 


Figure  42:  Comparison  of  results  for  G  computed  with  Discretizations  1  and  3. 
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Time  History,  Plane  Strain  Model 


Figure  43:  Time  history  for  Ki  and  Kn  using  coarse  mesh  and  Pc  =  0  +  1,  Pc  =  0  +  (2, 2, 1) 
(Discretizations  3  and  4). 


Time  History,  Plane  Strain  Model 


Figure  44:  Time  history  for  G  using  coarse  mesh  and  pc  =  0  + 1,  Pc  —  0  +  (2, 2, 1)  (Discretiza¬ 
tions  3  and  4). 
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Time  History,  Plane  Strain  Model 


Figure  45:  Time  history  for  G  using  coarse  mesh  and  Pc  =  0  + 1,  Pc  =  0  +  (2, 2, 1)  (Discretiza¬ 
tions  3  and  4).  The  tic  marks  of  the  x-axis  are  placed  exactly  when  the  crack  front  crossed 
the  boundary  between  two  elements. 

13.4  Edge  Cracked  Plate  Under  Asymmetric  Impact 

This  example  also  shows  the  excellent  agreement  between  solutions  obtained  with  PHLEX- 
crack  and  analytical  solutions.  The  problem  is  a  plate  containing  an  edge  crack  (Figure  46). 
Initially,  the  plate  is  stress-free  and  at  rest.  The  asymmetric  impact  of  a  projectile  on  the 
cracked  edge  is  simulated  as  a  normal  velocity  imposed  suddenly  at  t=0  on  the  side  of  the 
plate.  The  problem  and  its  analytical  solution  were  presented  by  Lee  and  Freund  [26].  Under 
the  current  loading,  the  crack  is  subjected  to  negative  mode  I  stress  intensity  factor  and, 
therefore,  will  be  stationary;  the  loading  is  predominantly  mode  II. 

The  problem  was  solved  using  both  the  finite  element  method  and  the  EPM  method. 
Figure  46  shows  both  models.  Note  that  in  the  finite  element  model,  successive  levels  of 
mesh  refinement  were  needed  at  the  crack  tip  in  order  to  obtain  accurate  results.  For  the 
EPM  model,  singular  functions  were  used  at  the  crack  tip.  Model  parameters  and  material 
properties  were  chosen  as  those  used  by  Organ  [44]: 

•  Young’s  modulus  =  E  =  200  x  lO^N/m^ 


•  Poisson’s  ratio  =  0.25 

•  Mass  density  =  7833  Kg/m^ 
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•  K  =  Kosolov  constant  =  3  —  Au  =  2  (for  plane  strain) 

•  Impact  speed  —  vq  =  16.5  m/sec 

•  Crack  length  =  oq  =  1.0  m 

•  Plate  dimensions  =  4ao  x  6ao 

Figure  47  presents  the  time  histories  of  mode  I  and  II  stress  intensity  factors  for  all  three 
solutions:  Analytical,  finite  element,  and  EPM.  The  stress  intensity  factors  were  normalized 
by 


-Evo^/oo/tt 

2(1  - 


-0.17939  X  10^ 


where 


K  +  l  IJl 
K  —  ly  p 


5535.3m/ s, 


is  the  dilatational  wave  speed.  Time  was  normalized  with 

ao/cd  =  18.066  x  10"® 

The  two  graphs  of  Figure  47  clearly  show  the  excellent  agreement  among  all  three 
solutions.  In  this  example,  the  natural  and  simplified  inclusion  of  singular  functions  with 
the  EPM  approximation  and  the  ability  to  place  the  crack  anywhere  in  the  model  without 
modifying  or  refining  it,  are  important  advantages  of  EPM  over  FE  method.  The  advantages 
of  the  EPM  method,  however,  will  become  more  evident  in  propagation  problems  where  the 
crack  can  freely  and  arbitrarily  advance  inside  the  domain  without  the  need  for  remeshing. 


13.5  Welded  T-Joint  with  a  Crack 

In  this  section  we  present  a  truly  3D  example  on  using  EPM  for  crack  propagation.  The 
example,  as  shown  in  Figure  48,  is  a  beam  with  welded  cross-section  (T-section).  An  initial 
half-penny  crack  is  placed  longitudinally  between  the  weld  and  the  web  as  shown  in  the  figure. 
The  crack  propagates  due  to  an  impact  couple  loading  applied  at  time=0  at  the  end  of  the 
beam.  Dynamic  waves  travel  through  the  body.  Once  it  reach  the  crack  area,  stresses  increase 
sharply  so  that  crack  is  propagated.  It  is  assumed  that  both  domain  and  loading  are  symmetric 
with  respect  to  the  other  end  of  the  beam  and,  therefore,  only  half  of  the  domain  is  analyzed. 
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(a)  EPM 


(b)  EPM  in  vicinity  of 
crack  tip 


L( 


(c)  FE 


(d)  FE  in  vicinity  of 
crack  tip 


Figure  46:  The  grids  used:  (a,b)  for  EPM,  and  (c,d)  for  FE  analyses. 


Time  History  of  SIF’s 


Normalized  Time 


(a)  Mode  I 


Time  History  of  SIF’s 


Normalized  Time 


(b)  Mode  II 

Figure  47:  Comparisons  of  Stress  Intensity  Factors  obtained  from  EPM,  FE,  and  analytical 
solutions. 


Figure  48:  Welded  T-section  beam  with  a  crack. 
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Both  the  web  and  the  flange  are  made  of  the  same  linear  elastic  material  with  the 
following  properties: 


•  Young’s  modulus  =  200  x  10® 

•  Poisson’s  ratio  =  0.3 

•  Mass  density  =  7833 


The  weld  material  is  assumed  to  be  ten  times  stiffer  than  the  latter  (i.  e.,  Young  s  modulus 
200  X  10^°).  The  Freund  propagation  model  is  used  to  advance  the  crack  (see  section  6).  The 
dynamic  fracture  properties  and  parameters  used  to  propagate  the  crack  are  as  follows: 


•  Least  squares  fit  extraction  domain:  Cylinder  of  radius  0.5  and  length  0.25 

•  Number  of  integration  points  in  the  r,  9  and  2:  directions:  n  =  (5, 10, 3),  respectively. 

•  Type  of  D  matrix  is  material  matrix. 

•  Power  for  the  weighting  function  =  3 

•  Dynamic  fracture  toughness  in  a  pure  mode  I,  Kjd  =  75000 

•  Crack  speed  limit  =  1212 


The  couple  is  applied  as  two  equal  and  opposite  impact  forces  at  the  two  corners  of  the  edge 
cross  section  with  a  value  of  5  x  10®  each. 

The  initial  crack  surface  definition  is  shown  in  Figure  48.  Note  that  the  web  is  not 
in  contact  with  the  flange,  but  rather  a  small  gap  exists  between  the  two.  The  domain 
is  discretized  for  EPM  as  shown  in  Figure  49.  The  figure  shows  the  integration  cells;  the 
meshless  nodes  are  located  at  the  corners  of  the  these  cells.  Note  that  the  grid  is  made  finer 
around  the  crack  area.  Linear  approximation  is  used  over  all  the  domain. 

A  transient  dynamic  analysis  is  performed  on  the  model  described  above.  Newmark 
method  is  used  to  march  the  solution  over  time.  Figures  50,  51,  and  52  represent  the  crack 
surface  at  times  0.0015,  0.0020,  and  0.0030  respectively. 

13.6  Welded  T-Joint  with  a  Fully  Embedded  Crack 

This  example  shows  identical  domain  as  in  the  previous  example,  but  the  initial  crack  has 
been  changed  to  full  disc  and  moved  to  the  inside  of  the  domain.  Fig.  53  shows  the  position 


80 


z 

Figure  49;  EPM  discretization  for  welded  T-section  beam. 

of  the  crack,  and  Fig.  54  shows  the  crack  at  the  time  0.0021.  This  crack  has  been  formed  due 
to  5  consecutive  shock  waves,  between  them  the  speed  of  crack  front  slowed  down  and  even 
stopped  (this  is  visible  as  darker  bands  of  edges  inside  the  crack  surface).  History  of  crack 
propagation  is  shown  in  Fig.  55. 

This  example  shows  also  start  of  numerical  instability  in  the  crack.  Details  of  this 
instability,  and  practical  means  to  control  it  have  been  discussed  in  Section.  12. 


Figure  51:  Crack  surface  of  welded  T-section  example  at  time  0.0020. 
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Figure  52:  Crack  surface  of  welded  T-section  example  at  time  0.0030. 
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Figure  54:  Top  view  of  the  penny  shape  crack  in  a  welded  T-section  beam.  (Shock  waves  are 
coming  along  Z-axis) 
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13.7  Section  of  Infinite  Cylinder 


This  example  shows  the  capability  of  the  code  to  perform  analysis  using  classical  finite  elements 
and  EPM  ICell  mesh.  The  initial  model  consists  of  half  of  the  cylinder,  with  a  rib  modeled 
using  shell  elements  (QUAD4).  Appropriate  symmetry  boundary  conditions  are  applied  to 
the  shell  and  EPM  elements.  The  domain  and  boundary  conditions  are  shown  in  Fig.  56. 


Figure  56:  Domain  and  boundary  conditions  for  the  infinite  cylinder  problem. 

The  physical  and  material  parameters  for  this  problem  are  meant  to  simulate  a  steel 
shell  but  are  non-dimensionalized. 

•  Geometric  properties: 


-  Shell  radius  =1.0 

-  Shell  and  rib  thickness  =  0.01 

-  Axial  width  =  0.175 
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—  Rib  inner  radius  =  0.825 

•  Material  properties: 

—  Young’s  modulus  £'  =  98.125 
—  Poisson’s  ratio  1/  =  0.3 

-  Mass  density  p  — 7.85 

—  Dynamic  fracture  toughness  Kjp  =  1.455  x  10  ® 

-  Rayleigh  wave  speed  Cr  =  2.0335 


The  load  is  applied  at  time  t  =  0  as  shown  in  Fig.  56.  The  point  loads  correspond  to 
a  uniform  pressure  of  magnitude  10”^  applied  to  the  faces  of  the  elements  shown  in  Fig.  56. 

The  following  sequence  of  problems  has  been  solved  for  this  model: 

1  Coarse  shell  mesh  (18  x  2  x  2  linear  elements,  i.e.  18  elements  along  the  arc,  2  in  radial 
and  2  in  the  Z  direction) 

2  Fine  shell  mesh  (54  x  6  x  6  linear  elements) 

3  Mixed  mesh  with  part  of  the  rib  modeled  using  EPM  approximation  (Figs.  57  and  58). 
The  shell  elements  are  linear  while  the  ICells  are  quadratic  in  the  X  and  Y  directions 
and  linear  in  the  Z  direction. 

4  Mixed  mesh  with  crack  (Fig.  59).  The  same  p-order  as  in  the  previous  problem  is  used 
here. 


The  dilatational  wave  speed  is  given  by 


Cd  = 


K  \  p 

K  —  1  p 


4.102 


where  k  =  3  —  Au  and  p  = 

In  problems  3  and  4,  the  time  step  was  chosen  equal  to  At  =  6  x  10^  which  is  about 
half  of  hminICd  where  hmin  is  the  smallest  element  size  in  the  mesh. 

In  the  case  of  problem  4  the  crack  speed  is  prescribed  as 

C{t)  =  Hit  -  0.06)^  =  H{t  -  0.06)0.3389 


By  comparing  results  of  shell-only  models  it  has  been  determined  that  the  model  is 
sufficient  to  obtain  qualitative  results,  but  further  refinement  (i.e.  finer  model)  will  probably 
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result  in  better  values  of  the  stresses  and  crack  behavior.  The  results  will  not  be  presented  here 
in  details,  only  following  figures  are  enclosed  to  illustrate  the  results.  Fig.  60  shows  stresses 
at  the  time  0.00009  in  the  model  without  the  crack,  Figs.  61  and  62  show  displacements  and 
stresses  in  the  presence  of  stationary  crack  at  the  time  0.0001.  Figures  63,  64  and  65  show 
the  crack  growth  at  the  times  t  =  0.072,  0.102,  0.138,  respectively. 


Figure  57:  Infinite  cylinder  model  using  classical  finite  elements  and  EPM  ICells. 
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Figure  60:  Infinite  cylinder:  Stresses  in  the  model  without  the  crack,  time 
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Figure  61:  Infinite  cylinder:  X  Displacements  in  the  model  with  the  crack,  time 
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Figure  62:  Infinite  cylinder:  Stresses  in  the  model  with  the  crack,  time  =  0.0001 


Figure  63:  Infinite  cylinder:  Crack  growth  at  the  time  t  =  0.072. 
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Figure  64:  Infinite  cylinder:  Crack  growth  at  the  time  t  =  0.102. 


Figure  65:  Infinite  cylinder:  Crack  growth  at  the  time  t  =  0.138. 
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